Cantera  2.0
vcs_solve.cpp
1 /*!
2  * @file vcs_solve.h
3  * Header file for the internal class that holds 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 
11 
13 #include "vcs_Exception.h"
15 #include "cantera/equil/vcs_prob.h"
16 
18 #include "vcs_SpeciesProperties.h"
19 #include "vcs_species_thermo.h"
20 
21 #include "cantera/base/clockWC.h"
22 
23 #include <string>
24 #include <cstdlib>
25 #include <math.h>
26 
27 using namespace std;
28 
29 namespace VCSnonideal
30 {
31 
32 int vcs_timing_print_lvl = 1;
33 
34 VCS_SOLVE::VCS_SOLVE() :
35  NSPECIES0(0),
36  NPHASE0(0),
37  m_numSpeciesTot(0),
38  m_numElemConstraints(0),
39  m_numComponents(0),
40  m_numRxnTot(0),
41  m_numSpeciesRdc(0),
42  m_numRxnRdc(0),
43  m_numRxnMinorZeroed(0),
44  m_numPhases(0),
45  m_doEstimateEquil(0),
46  m_totalMolNum(0.0),
47  m_temperature(0.0),
48  m_pressurePA(0.0),
49  m_tolmaj(0.0),
50  m_tolmin(0.0),
51  m_tolmaj2(0.0),
52  m_tolmin2(0.0),
53  m_unitsState(VCS_DIMENSIONAL_G),
54  m_totalMoleScale(1.0),
55  m_useActCoeffJac(0),
56  m_totalVol(0.0),
57  m_Faraday_dim(Cantera::ElectronCharge * Cantera::Avogadro),
58  m_VCount(0),
59  m_debug_print_lvl(0),
60  m_timing_print_lvl(1),
61  m_VCS_UnitsFormat(VCS_UNITS_UNITLESS)
62 {
63 }
64 
65 // Initialize the sizes within the VCS_SOLVE object
66 
67 /*
68  * This resizes all of the internal arrays within the object. This routine
69  * operates in two modes. If all of the parameters are the same as it
70  * currently exists in the object, nothing is done by this routine; a quick
71  * exit is carried out and all of the data in the object persists.
72  *
73  * IF any of the parameters are different than currently exists in the
74  * object, then all of the data in the object must be redone. It may not
75  * be zeroed, but it must be redone.
76  *
77  * @param nspecies0 Number of species within the object
78  * @param nelements Number of element constraints within the problem
79  * @param nphase0 Number of phases defined within the problem.
80  *
81  */
82 void VCS_SOLVE::vcs_initSizes(const size_t nspecies0, const size_t nelements,
83  const size_t nphase0)
84 {
85 
86  if (NSPECIES0 != 0) {
87  if ((nspecies0 != NSPECIES0) || (nelements != m_numElemConstraints) || (nphase0 != NPHASE0)) {
89  } else {
90  return;
91  }
92  }
93 
94  NSPECIES0 = nspecies0;
95  NPHASE0 = nphase0;
96  m_numSpeciesTot = nspecies0;
97  m_numElemConstraints = nelements;
98  m_numComponents = nelements;
99 
100  string ser = "VCS_SOLVE: ERROR:\n\t";
101  if (nspecies0 <= 0) {
102  plogf("%s Number of species is nonpositive\n", ser.c_str());
103  throw vcsError("VCS_SOLVE()",
104  ser + " Number of species is nonpositive\n",
105  VCS_PUB_BAD);
106  }
107  if (nelements <= 0) {
108  plogf("%s Number of elements is nonpositive\n", ser.c_str());
109  throw vcsError("VCS_SOLVE()",
110  ser + " Number of species is nonpositive\n",
111  VCS_PUB_BAD);
112  }
113  if (nphase0 <= 0) {
114  plogf("%s Number of phases is nonpositive\n", ser.c_str());
115  throw vcsError("VCS_SOLVE()",
116  ser + " Number of species is nonpositive\n",
117  VCS_PUB_BAD);
118  }
119 
120  //vcs_priv_init(this);
121  m_VCS_UnitsFormat = VCS_UNITS_UNITLESS;
122 
123  /*
124  * We will initialize sc[] to note the fact that it needs to be
125  * filled with meaningful information.
126  */
127  m_stoichCoeffRxnMatrix.resize(nspecies0, nelements, 0.0);
128 
129  m_scSize.resize(nspecies0, 0.0);
130  m_spSize.resize(nspecies0, 1.0);
131 
132  m_SSfeSpecies.resize(nspecies0, 0.0);
133  m_feSpecies_new.resize(nspecies0, 0.0);
134  m_molNumSpecies_old.resize(nspecies0, 0.0);
135 
137 
138  m_deltaMolNumPhase.resize(nspecies0, nphase0, 0.0);
139  m_phaseParticipation.resize(nspecies0, nphase0, 0);
140  m_phasePhi.resize(nphase0, 0.0);
141 
142  m_molNumSpecies_new.resize(nspecies0, 0.0);
143 
144  m_deltaGRxn_new.resize(nspecies0, 0.0);
145  m_deltaGRxn_old.resize(nspecies0, 0.0);
146  m_deltaGRxn_Deficient.resize(nspecies0, 0.0);
147  m_deltaGRxn_tmp.resize(nspecies0, 0.0);
148  m_deltaMolNumSpecies.resize(nspecies0, 0.0);
149 
150  m_feSpecies_old.resize(nspecies0, 0.0);
151  m_elemAbundances.resize(nelements, 0.0);
152  m_elemAbundancesGoal.resize(nelements, 0.0);
153 
154  m_tPhaseMoles_old.resize(nphase0, 0.0);
155  m_tPhaseMoles_new.resize(nphase0, 0.0);
156  m_deltaPhaseMoles.resize(nphase0, 0.0);
157  m_TmpPhase.resize(nphase0, 0.0);
158  m_TmpPhase2.resize(nphase0, 0.0);
159 
160  m_formulaMatrix.resize(nelements, nspecies0);
161 
162  TPhInertMoles.resize(nphase0, 0.0);
163 
164  /*
165  * ind[] is an index variable that keep track of solution vector
166  * rotations.
167  */
168  m_speciesMapIndex.resize(nspecies0, 0);
169  m_speciesLocalPhaseIndex.resize(nspecies0, 0);
170  /*
171  * IndEl[] is an index variable that keep track of element vector
172  * rotations.
173  */
174  m_elementMapIndex.resize(nelements, 0);
175 
176  /*
177  * ir[] is an index vector that keeps track of the irxn to species
178  * mapping. We can't fill it in until we know the number of c
179  * components in the problem
180  */
181  m_indexRxnToSpecies.resize(nspecies0, 0);
182 
183  /* Initialize all species to be major species */
184  //m_rxnStatus.resize(nspecies0, 1);
185  m_speciesStatus.resize(nspecies0, 1);
186 
187  m_SSPhase.resize(2*nspecies0, 0);
188  m_phaseID.resize(nspecies0, 0);
189 
190  m_numElemConstraints = nelements;
191 
192  m_elementName.resize(nelements, std::string(""));
193  m_speciesName.resize(nspecies0, std::string(""));
194 
195  m_elType.resize(nelements, VCS_ELEM_TYPE_ABSPOS);
196 
197  m_elementActive.resize(nelements, 1);
198  /*
199  * Malloc space for activity coefficients for all species
200  * -> Set it equal to one.
201  */
202  m_actConventionSpecies.resize(nspecies0, 0);
203  m_phaseActConvention.resize(nphase0, 0);
204  m_lnMnaughtSpecies.resize(nspecies0, 0.0);
205  m_actCoeffSpecies_new.resize(nspecies0, 1.0);
206  m_actCoeffSpecies_old.resize(nspecies0, 1.0);
207  m_wtSpecies.resize(nspecies0, 0.0);
208  m_chargeSpecies.resize(nspecies0, 0.0);
209  m_speciesThermoList.resize(nspecies0, (VCS_SPECIES_THERMO*)0);
210 
211  /*
212  * Malloc Phase Info
213  */
214  m_VolPhaseList.resize(nphase0, 0);
215  for (size_t iph = 0; iph < nphase0; iph++) {
216  m_VolPhaseList[iph] = new vcs_VolPhase(this);
217  }
218 
219  /*
220  * For Future expansion
221  */
222  m_useActCoeffJac = true;
223  if (m_useActCoeffJac) {
224  m_dLnActCoeffdMolNum.resize(nspecies0, nspecies0, 0.0);
225  }
226 
227  m_PMVolumeSpecies.resize(nspecies0, 0.0);
228 
229  /*
230  * Malloc space for counters kept within vcs
231  *
232  */
233  m_VCount = new VCS_COUNTERS();
235 
236  if (vcs_timing_print_lvl == 0) {
237  m_timing_print_lvl = 0;
238  }
239 
240  return;
241 
242 }
243 /****************************************************************************/
244 
245 // Destructor
247 {
249 }
250 /*****************************************************************************/
251 
252 // Delete memory that isn't just resizeable STL containers
253 /*
254  * This gets called by the destructor or by InitSizes().
255  */
257 {
258  size_t j;
259  size_t nspecies = m_numSpeciesTot;
260 
261  for (j = 0; j < m_numPhases; j++) {
262  delete m_VolPhaseList[j];
263  m_VolPhaseList[j] = 0;
264  }
265 
266  for (j = 0; j < nspecies; j++) {
267  delete m_speciesThermoList[j];
268  m_speciesThermoList[j] = 0;
269  }
270 
271  delete m_VCount;
272  m_VCount = 0;
273 
274  NSPECIES0 = 0;
275  NPHASE0 = 0;
277  m_numComponents = 0;
278 }
279 /*****************************************************************************/
280 
281 // Solve an equilibrium problem
282 /*
283  * This is the main interface routine to the equilibrium solver
284  *
285  * Input:
286  * @param vprob Object containing the equilibrium Problem statement
287  *
288  * @param ifunc Determines the operation to be done: Valid values:
289  * 0 -> Solve a new problem by initializing structures
290  * first. An initial estimate may or may not have
291  * been already determined. This is indicated in the
292  * VCS_PROB structure.
293  * 1 -> The problem has already been initialized and
294  * set up. We call this routine to resolve it
295  * using the problem statement and
296  * solution estimate contained in
297  * the VCS_PROB structure.
298  * 2 -> Don't solve a problem. Destroy all the private
299  * structures.
300  *
301  * @param ipr Printing of results
302  * ipr = 1 -> Print problem statement and final results to
303  * standard output
304  * 0 -> don't report on anything
305  * @param ip1 Printing of intermediate results
306  * ip1 = 1 -> Print intermediate results.
307  * = 0 -> No intermediate results printing
308  *
309  * @param maxit Maximum number of iterations for the algorithm
310  *
311  * Output:
312  *
313  * @return
314  * nonzero value: failure to solve the problem at hand.
315  * zero : success
316  */
317 int VCS_SOLVE::vcs(VCS_PROB* vprob, int ifunc, int ipr, int ip1, int maxit)
318 {
319  int retn = 0, iconv = 0;
320  size_t nspecies0, nelements0, nphase0;
321  Cantera::clockWC tickTock;
322 
323  int iprintTime = std::max(ipr, ip1);
324  if (m_timing_print_lvl < iprintTime) {
325  iprintTime = m_timing_print_lvl ;
326  }
327 
328  if (ifunc > 2) {
329  plogf("vcs: Unrecognized value of ifunc, %d: bailing!\n",
330  ifunc);
331  return VCS_PUB_BAD;
332  }
333 
334  if (ifunc == 0) {
335  /*
336  * This function is called to create the private data
337  * using the public data.
338  */
339  nspecies0 = vprob->nspecies + 10;
340  nelements0 = vprob->ne;
341  nphase0 = vprob->NPhase;
342 
343  vcs_initSizes(nspecies0, nelements0, nphase0);
344 
345  if (retn != 0) {
346  plogf("vcs_priv_alloc returned a bad status, %d: bailing!\n",
347  retn);
348  return retn;
349  }
350  /*
351  * This function is called to copy the public data
352  * and the current problem specification
353  * into the current object's data structure.
354  */
355  retn = vcs_prob_specifyFully(vprob);
356  if (retn != 0) {
357  plogf("vcs_pub_to_priv returned a bad status, %d: bailing!\n",
358  retn);
359  return retn;
360  }
361  /*
362  * Prep the problem data
363  * - adjust the identity of any phases
364  * - determine the number of components in the problem
365  */
366  retn = vcs_prep_oneTime(ip1);
367  if (retn != 0) {
368  plogf("vcs_prep_oneTime returned a bad status, %d: bailing!\n",
369  retn);
370  return retn;
371  }
372  }
373  if (ifunc == 1) {
374  /*
375  * This function is called to copy the current problem
376  * into the current object's data structure.
377  */
378  retn = vcs_prob_specify(vprob);
379  if (retn != 0) {
380  plogf("vcs_prob_specify returned a bad status, %d: bailing!\n",
381  retn);
382  return retn;
383  }
384  }
385  if (ifunc != 2) {
386  /*
387  * Prep the problem data for this particular instantiation of
388  * the problem
389  */
390  retn = vcs_prep();
391  if (retn != VCS_SUCCESS) {
392  plogf("vcs_prep returned a bad status, %d: bailing!\n", retn);
393  return retn;
394  }
395 
396  /*
397  * Check to see if the current problem is well posed.
398  */
399  if (!vcs_wellPosed(vprob)) {
400  plogf("vcs has determined the problem is not well posed: Bailing\n");
401  return VCS_PUB_BAD;
402  }
403 
404  /*
405  * Once we have defined the global internal data structure defining
406  * the problem, then we go ahead and solve the problem.
407  *
408  * (right now, all we do is solve fixed T, P problems.
409  * Methods for other problem types will go in at this level.
410  * For example, solving for fixed T, V problems will involve
411  * a 2x2 Newton's method, using loops over vcs_TP() to
412  * calculate the residual and Jacobian)
413  */
414  iconv = vcs_TP(ipr, ip1, maxit, vprob->T, vprob->PresPA);
415 
416  /*
417  * If requested to print anything out, go ahead and do so;
418  */
419  if (ipr > 0) {
420  vcs_report(iconv);
421  }
422  /*
423  * Copy the results of the run back to the VCS_PROB structure,
424  * which is returned to the user.
425  */
426  vcs_prob_update(vprob);
427  }
428 
429  /*
430  * Report on the time if requested to do so
431  */
432  double te = tickTock.secondsWC();
433  m_VCount->T_Time_vcs += te;
434  if (iprintTime > 0) {
436  }
437  /*
438  * Now, destroy the private data, if requested to do so
439  *
440  * FILL IN
441  */
442 
443  if (iconv < 0) {
444  plogf("ERROR: FAILURE its = %d!\n", m_VCount->Its);
445  } else if (iconv == 1) {
446  plogf("WARNING: RANGE SPACE ERROR encountered\n");
447  }
448  return iconv;
449 }
450 /*****************************************************************************/
451 
452 // Fully specify the problem to be solved using VCS_PROB
453 /*
454  * Use the contents of the VCS_PROB to specify the contents of the
455  * private data, VCS_SOLVE.
456  *
457  * @param pub Pointer to VCS_PROB that will be used to
458  * initialize the current equilibrium problem
459  */
461 {
462  vcs_VolPhase* Vphase = 0;
463  const char* ser =
464  "vcs_pub_to_priv ERROR :ill defined interface -> bailout:\n\t";
465 
466  /*
467  * First Check to see whether we have room for the current problem
468  * size
469  */
470  size_t nspecies = pub->nspecies;
471  if (NSPECIES0 < nspecies) {
472  plogf("%sPrivate Data is dimensioned too small\n", ser);
473  return VCS_PUB_BAD;
474  }
475  size_t nph = pub->NPhase;
476  if (NPHASE0 < nph) {
477  plogf("%sPrivate Data is dimensioned too small\n", ser);
478  return VCS_PUB_BAD;
479  }
480  size_t nelements = pub->ne;
481  if (m_numElemConstraints < nelements) {
482  plogf("%sPrivate Data is dimensioned too small\n", ser);
483  return VCS_PUB_BAD;
484  }
485 
486  /*
487  * OK, We have room. Now, transfer the integer numbers
488  */
489  m_numElemConstraints = nelements;
490  m_numSpeciesTot = nspecies;
492  /*
493  * nc = number of components -> will be determined later.
494  * but set it to its maximum possible value here.
495  */
496  m_numComponents = nelements;
497  /*
498  * m_numRxnTot = number of noncomponents, also equal to the
499  * number of reactions
500  * Note, it's possible that the number of elements is greater than
501  * the number of species. In that case set the number of reactions
502  * to zero.
503  */
504  if (nelements > nspecies) {
505  m_numRxnTot = 0;
506  } else {
507  m_numRxnTot = nspecies - nelements;
508  }
510  /*
511  * number of minor species rxn -> all species rxn are major at the start.
512  */
514  /*
515  * NPhase = number of phases
516  */
517  m_numPhases = nph;
518 
519 #ifdef DEBUG_MODE
521 #else
523 #endif
524 
525  /*
526  * FormulaMatrix[] -> Copy the formula matrix over
527  */
528  for (size_t i = 0; i < nspecies; i++) {
529  bool nonzero = false;
530  for (size_t j = 0; j < nelements; j++) {
531  if (pub->FormulaMatrix[j][i] != 0.0) {
532  nonzero = true;
533  }
534  m_formulaMatrix[j][i] = pub->FormulaMatrix[j][i];
535  }
536  if (!nonzero) {
537  plogf("vcs_prob_specifyFully:: species %d %s has a zero formula matrix!\n", i,
538  pub->SpName[i].c_str());
539  return VCS_PUB_BAD;
540  }
541  }
542 
543  /*
544  * Copy over the species molecular weights
545  */
546  vcs_vdcopy(m_wtSpecies, pub->WtSpecies, nspecies);
547 
548  /*
549  * Copy over the charges
550  */
551  vcs_vdcopy(m_chargeSpecies, pub->Charge, nspecies);
552 
553  /*
554  * Malloc and Copy the VCS_SPECIES_THERMO structures
555  *
556  */
557  for (size_t kspec = 0; kspec < nspecies; kspec++) {
558  if (m_speciesThermoList[kspec] != NULL) {
559  delete m_speciesThermoList[kspec];
560  }
561  VCS_SPECIES_THERMO* spf = pub->SpeciesThermo[kspec];
562  m_speciesThermoList[kspec] = spf->duplMyselfAsVCS_SPECIES_THERMO();
563  if (m_speciesThermoList[kspec] == NULL) {
564  plogf(" duplMyselfAsVCS_SPECIES_THERMO returned an error!\n");
565  return VCS_PUB_BAD;
566  }
567  }
568 
569  /*
570  * Copy the species unknown type
571  */
573  VCS_DATA_PTR(pub->SpeciesUnknownType), nspecies);
574 
575  /*
576  * iest => Do we have an initial estimate of the species mole numbers ?
577  */
578  m_doEstimateEquil = pub->iest;
579 
580  /*
581  * w[] -> Copy the equilibrium mole number estimate if it exists.
582  */
583  if (pub->w.size() != 0) {
584  vcs_vdcopy(m_molNumSpecies_old, pub->w, nspecies);
585  } else {
586  m_doEstimateEquil = -1;
587  vcs_dzero(VCS_DATA_PTR(m_molNumSpecies_old), nspecies);
588  }
589 
590  /*
591  * Formulate the Goal Element Abundance Vector
592  */
593  if (pub->gai.size() != 0) {
594  for (size_t i = 0; i < nelements; i++) {
595  m_elemAbundancesGoal[i] = pub->gai[i];
596  if (pub->m_elType[i] == VCS_ELEM_TYPE_LATTICERATIO) {
597  if (m_elemAbundancesGoal[i] < 1.0E-10) {
598  m_elemAbundancesGoal[i] = 0.0;
599  }
600  }
601  }
602  } else {
603  if (m_doEstimateEquil == 0) {
604  double sum = 0;
605  for (size_t j = 0; j < nelements; j++) {
606  m_elemAbundancesGoal[j] = 0.0;
607  for (size_t kspec = 0; kspec < nspecies; kspec++) {
609  sum += m_molNumSpecies_old[kspec];
611  }
612  }
613  if (pub->m_elType[j] == VCS_ELEM_TYPE_LATTICERATIO) {
614  if (m_elemAbundancesGoal[j] < 1.0E-10 * sum) {
615  m_elemAbundancesGoal[j] = 0.0;
616  }
617  }
618  }
619  } else {
620  plogf("%sElement Abundances, m_elemAbundancesGoal[], not specified\n", ser);
621  return VCS_PUB_BAD;
622  }
623  }
624 
625  /*
626  * zero out values that will be filled in later
627  */
628  /*
629  * TPhMoles[] -> Untouched here. These will be filled in vcs_prep.c
630  * TPhMoles1[]
631  * DelTPhMoles[]
632  *
633  *
634  * T, Pres, copy over here
635  */
636  if (pub->T > 0.0) {
637  m_temperature = pub->T;
638  } else {
639  m_temperature = 293.15;
640  }
641  if (pub->PresPA > 0.0) {
642  m_pressurePA = pub->PresPA;
643  } else {
645  }
646  /*
647  * TPhInertMoles[] -> must be copied over here
648  */
649  for (size_t iph = 0; iph < nph; iph++) {
650  Vphase = pub->VPhaseList[iph];
651  TPhInertMoles[iph] = Vphase->totalMolesInert();
652  }
653 
654  /*
655  * if__ : Copy over the units for the chemical potential
656  */
658 
659  /*
660  * tolerance requirements -> copy them over here and later
661  */
662  m_tolmaj = pub->tolmaj;
663  m_tolmin = pub->tolmin;
664  m_tolmaj2 = 0.01 * m_tolmaj;
665  m_tolmin2 = 0.01 * m_tolmin;
666 
667  /*
668  * m_speciesIndexVector[] is an index variable that keep track
669  * of solution vector rotations.
670  */
671  for (size_t i = 0; i < nspecies; i++) {
672  m_speciesMapIndex[i] = i;
673  }
674 
675  /*
676  * IndEl[] is an index variable that keep track of element vector
677  * rotations.
678  */
679  for (size_t i = 0; i < nelements; i++) {
680  m_elementMapIndex[i] = i;
681  }
682 
683  /*
684  * Define all species to be major species, initially.
685  */
686  for (size_t i = 0; i < nspecies; i++) {
687  // m_rxnStatus[i] = VCS_SPECIES_MAJOR;
689  }
690  /*
691  * PhaseID: Fill in the species to phase mapping
692  * -> Check for bad values at the same time.
693  */
694  if (pub->PhaseID.size() != 0) {
695  std::vector<size_t> numPhSp(nph, 0);
696  for (size_t kspec = 0; kspec < nspecies; kspec++) {
697  size_t iph = pub->PhaseID[kspec];
698  if (iph >= nph) {
699  plogf("%sSpecies to Phase Mapping, PhaseID, has a bad value\n",
700  ser);
701  plogf("\tPhaseID[%d] = %d\n", kspec, iph);
702  plogf("\tAllowed values: 0 to %d\n", nph - 1);
703  return VCS_PUB_BAD;
704  }
705  m_phaseID[kspec] = pub->PhaseID[kspec];
706  m_speciesLocalPhaseIndex[kspec] = numPhSp[iph];
707  numPhSp[iph]++;
708  }
709  for (size_t iph = 0; iph < nph; iph++) {
710  Vphase = pub->VPhaseList[iph];
711  if (numPhSp[iph] != Vphase->nSpecies()) {
712  plogf("%sNumber of species in phase %d, %s, doesn't match\n",
713  ser, iph, Vphase->PhaseName.c_str());
714  return VCS_PUB_BAD;
715  }
716  }
717  } else {
718  if (m_numPhases == 1) {
719  for (size_t kspec = 0; kspec < nspecies; kspec++) {
720  m_phaseID[kspec] = 0;
721  m_speciesLocalPhaseIndex[kspec] = kspec;
722  }
723  } else {
724  plogf("%sSpecies to Phase Mapping, PhaseID, is not defined\n", ser);
725  return VCS_PUB_BAD;
726  }
727  }
728 
729  /*
730  * Copy over the element types
731  */
732  m_elType.resize(nelements, VCS_ELEM_TYPE_ABSPOS);
733  m_elementActive.resize(nelements, 1);
734 
735  /*
736  * Copy over the element names and types
737  */
738  for (size_t i = 0; i < nelements; i++) {
739  m_elementName[i] = pub->ElName[i];
740  m_elType[i] = pub->m_elType[i];
741  m_elementActive[i] = pub->ElActive[i];
742  if (!strncmp(m_elementName[i].c_str(), "cn_", 3)) {
744  if (pub->m_elType[i] != VCS_ELEM_TYPE_CHARGENEUTRALITY) {
745  plogf("we have an inconsistency!\n");
746  exit(EXIT_FAILURE);
747  }
748  }
749  }
750 
751  for (size_t i = 0; i < nelements; i++) {
753  if (m_elemAbundancesGoal[i] != 0.0) {
754  if (fabs(m_elemAbundancesGoal[i]) > 1.0E-9) {
755  plogf("Charge neutrality condition %s is signicantly nonzero, %g. Giving up\n",
756  m_elementName[i].c_str(), m_elemAbundancesGoal[i]);
757  exit(EXIT_FAILURE);
758  } else {
759  if (m_debug_print_lvl >= 2) {
760  plogf("Charge neutrality condition %s not zero, %g. Setting it zero\n",
761  m_elementName[i].c_str(), m_elemAbundancesGoal[i]);
762  }
763  m_elemAbundancesGoal[i] = 0.0;
764  }
765 
766  }
767  }
768  }
769 
770  /*
771  * Copy over the species names
772  */
773  for (size_t i = 0; i < nspecies; i++) {
774  m_speciesName[i] = pub->SpName[i];
775  }
776  /*
777  * Copy over all of the phase information
778  * Use the object's assignment operator
779  */
780  for (size_t iph = 0; iph < nph; iph++) {
781  *(m_VolPhaseList[iph]) = *(pub->VPhaseList[iph]);
782  /*
783  * Fix up the species thermo pointer in the vcs_SpeciesThermo object
784  * It should point to the species thermo pointer in the private
785  * data space.
786  */
787  Vphase = m_VolPhaseList[iph];
788  for (size_t k = 0; k < Vphase->nSpecies(); k++) {
789  vcs_SpeciesProperties* sProp = Vphase->speciesProperty(k);
790  size_t kT = Vphase->spGlobalIndexVCS(k);
791  sProp->SpeciesThermo = m_speciesThermoList[kT];
792  }
793  }
794 
795  /*
796  * Specify the Activity Convention information
797  */
798  for (size_t iph = 0; iph < nph; iph++) {
799  Vphase = m_VolPhaseList[iph];
801  if (Vphase->p_activityConvention != 0) {
802  /*
803  * We assume here that species 0 is the solvent.
804  * The solvent isn't on a unity activity basis
805  * The activity for the solvent assumes that the
806  * it goes to one as the species mole fraction goes to
807  * one; i.e., it's really on a molarity framework.
808  * So SpecLnMnaught[iSolvent] = 0.0, and the
809  * loop below starts at 1, not 0.
810  */
811  size_t iSolvent = Vphase->spGlobalIndexVCS(0);
812  double mnaught = m_wtSpecies[iSolvent] / 1000.;
813  for (size_t k = 1; k < Vphase->nSpecies(); k++) {
814  size_t kspec = Vphase->spGlobalIndexVCS(k);
816  m_lnMnaughtSpecies[kspec] = log(mnaught);
817  }
818  }
819  }
820 
821  /*
822  * Copy the title info
823  */
824  if (pub->Title.size() == 0) {
825  m_title = "Unspecified Problem Title";
826  } else {
827  m_title = pub->Title;
828  }
829 
830  /*
831  * Copy the volume info
832  */
833  m_totalVol = pub->Vol;
834  if (m_PMVolumeSpecies.size() != 0) {
835  vcs_dcopy(VCS_DATA_PTR(m_PMVolumeSpecies), VCS_DATA_PTR(pub->VolPM), nspecies);
836  }
837 
838  /*
839  * Return the success flag
840  */
841  return VCS_SUCCESS;
842 }
843 /*****************************************************************************/
844 
845 // Specify the problem to be solved using VCS_PROB, incrementally
846 /*
847  * Use the contents of the VCS_PROB to specify the contents of the
848  * private data, VCS_SOLVE.
849  *
850  * It's assumed we are solving the same problem.
851  *
852  * @param pub Pointer to VCS_PROdB that will be used to
853  * initialize the current equilibrium problem
854  */
856 {
857  size_t kspec, k, i, j, iph;
858  string yo("vcs_prob_specify ERROR: ");
859  int retn = VCS_SUCCESS;
860  bool status_change = false;
861 
862  m_temperature = pub->T;
863  m_pressurePA = pub->PresPA;
865  m_doEstimateEquil = pub->iest;
866 
867  m_totalVol = pub->Vol;
868 
869  m_tolmaj = pub->tolmaj;
870  m_tolmin = pub->tolmin;
871  m_tolmaj2 = 0.01 * m_tolmaj;
872  m_tolmin2 = 0.01 * m_tolmin;
873 
874  for (kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
875  k = m_speciesMapIndex[kspec];
876  m_molNumSpecies_old[kspec] = pub->w[k];
877  m_molNumSpecies_new[kspec] = pub->mf[k];
878  m_feSpecies_old[kspec] = pub->m_gibbsSpecies[k];
879  }
880 
881  /*
882  * Transfer the element abundance goals to the solve object
883  */
884  for (i = 0; i < m_numElemConstraints; i++) {
885  j = m_elementMapIndex[i];
886  m_elemAbundancesGoal[i] = pub->gai[j];
887  }
888 
889  /*
890  * Try to do the best job at guessing at the title
891  */
892  if (pub->Title.size() == 0) {
893  if (m_title.size() == 0) {
894  m_title = "Unspecified Problem Title";
895  }
896  } else {
897  m_title = pub->Title;
898  }
899 
900  /*
901  * Copy over the phase information.
902  * -> For each entry in the phase structure, determine
903  * if that entry can change from its initial value
904  * Either copy over the new value or create an error
905  * condition.
906  */
907 
908  for (iph = 0; iph < m_numPhases; iph++) {
909  vcs_VolPhase* vPhase = m_VolPhaseList[iph];
910  vcs_VolPhase* pub_phase_ptr = pub->VPhaseList[iph];
911 
912  if (vPhase->VP_ID_ != pub_phase_ptr->VP_ID_) {
913  plogf("%sPhase numbers have changed:%d %d\n", yo.c_str(),
914  vPhase->VP_ID_, pub_phase_ptr->VP_ID_);
915  retn = VCS_PUB_BAD;
916  }
917 
918  if (vPhase->m_singleSpecies != pub_phase_ptr->m_singleSpecies) {
919  plogf("%sSingleSpecies value have changed:%d %d\n", yo.c_str(),
920  vPhase->m_singleSpecies,
921  pub_phase_ptr->m_singleSpecies);
922  retn = VCS_PUB_BAD;
923  }
924 
925  if (vPhase->m_gasPhase != pub_phase_ptr->m_gasPhase) {
926  plogf("%sGasPhase value have changed:%d %d\n", yo.c_str(),
927  vPhase->m_gasPhase,
928  pub_phase_ptr->m_gasPhase);
929  retn = VCS_PUB_BAD;
930  }
931 
932  vPhase->m_eqnState = pub_phase_ptr->m_eqnState;
933 
934  if (vPhase->nSpecies() != pub_phase_ptr->nSpecies()) {
935  plogf("%sNVolSpecies value have changed:%d %d\n", yo.c_str(),
936  vPhase->nSpecies(),
937  pub_phase_ptr->nSpecies());
938  retn = VCS_PUB_BAD;
939  }
940 
941  if (vPhase->PhaseName != pub_phase_ptr->PhaseName) {
942  plogf("%sPhaseName value have changed:%s %s\n", yo.c_str(),
943  vPhase->PhaseName.c_str(),
944  pub_phase_ptr->PhaseName.c_str());
945  retn = VCS_PUB_BAD;
946  }
947 
948  if (vPhase->totalMolesInert() != pub_phase_ptr->totalMolesInert()) {
949  status_change = true;
950  }
951  /*
952  * Copy over the number of inert moles if it has changed.
953  */
954  TPhInertMoles[iph] = pub_phase_ptr->totalMolesInert();
955  vPhase->setTotalMolesInert(pub_phase_ptr->totalMolesInert());
956  if (TPhInertMoles[iph] > 0.0) {
957  vPhase->setExistence(2);
958  vPhase->m_singleSpecies = false;
959  }
960 
961  /*
962  * Copy over the interfacial potential
963  */
964  double phi = pub_phase_ptr->electricPotential();
965  vPhase->setElectricPotential(phi);
966  }
967 
968 
969  if (status_change) {
970  vcs_SSPhase();
971  }
972  /*
973  * Calculate the total number of moles in all phases.
974  */
975  vcs_tmoles();
976 
977  return retn;
978 }
979 /*****************************************************************************/
980 
981 // Transfer the results of the equilibrium calculation back to VCS_PROB
982 /*
983  * The VCS_PUB structure is returned to the user.
984  *
985  * @param pub Pointer to VCS_PROB that will get the results of the
986  * equilibrium calculation transfered to it.
987  */
989 {
990  size_t k1 = 0;
991 
992  vcs_tmoles();
995 
996  for (size_t i = 0; i < m_numSpeciesTot; ++i) {
997  /*
998  * Find the index of I in the index vector, m_speciesIndexVector[].
999  * Call it K1 and continue.
1000  */
1001  for (size_t j = 0; j < m_numSpeciesTot; ++j) {
1002  k1 = j;
1003  if (m_speciesMapIndex[j] == i) {
1004  break;
1005  }
1006  }
1007  /*
1008  * - Switch the species data back from K1 into I
1009  */
1011  pub->w[i] = m_molNumSpecies_old[k1];
1012  } else {
1013  pub->w[i] = 0.0;
1014  // plogf("voltage species = %g\n", m_molNumSpecies_old[k1]);
1015  }
1016  //pub->mf[i] = m_molNumSpecies_new[k1];
1017  pub->m_gibbsSpecies[i] = m_feSpecies_old[k1];
1018  pub->VolPM[i] = m_PMVolumeSpecies[k1];
1019  }
1020 
1021  pub->T = m_temperature;
1022  pub->PresPA = m_pressurePA;
1023  pub->Vol = m_totalVol;
1024  size_t kT = 0;
1025  for (size_t iph = 0; iph < pub->NPhase; iph++) {
1026  vcs_VolPhase* pubPhase = pub->VPhaseList[iph];
1027  vcs_VolPhase* vPhase = m_VolPhaseList[iph];
1028  pubPhase->setTotalMolesInert(vPhase->totalMolesInert());
1029  pubPhase->setTotalMoles(vPhase->totalMoles());
1030  pubPhase->setElectricPotential(vPhase->electricPotential());
1031  double sumMoles = pubPhase->totalMolesInert();
1032  pubPhase->setMoleFractionsState(vPhase->totalMoles(),
1033  VCS_DATA_PTR(vPhase->moleFractions()),
1035  const std::vector<double> & mfVector = pubPhase->moleFractions();
1036  for (size_t k = 0; k < pubPhase->nSpecies(); k++) {
1037  kT = pubPhase->spGlobalIndexVCS(k);
1038  pub->mf[kT] = mfVector[k];
1039  if (pubPhase->phiVarIndex() == k) {
1040  k1 = vPhase->spGlobalIndexVCS(k);
1041  double tmp = m_molNumSpecies_old[k1];
1042  if (! vcs_doubleEqual(pubPhase->electricPotential() , tmp)) {
1043  plogf("We have an inconsistency in voltage, %g, %g\n",
1044  pubPhase->electricPotential(), tmp);
1045  exit(EXIT_FAILURE);
1046  }
1047  }
1048 
1049  if (! vcs_doubleEqual(pub->mf[kT], vPhase->molefraction(k))) {
1050  plogf("We have an inconsistency in mole fraction, %g, %g\n",
1051  pub->mf[kT], vPhase->molefraction(k));
1052  exit(EXIT_FAILURE);
1053  }
1055  sumMoles += pub->w[kT];
1056  }
1057  }
1058  if (! vcs_doubleEqual(sumMoles, vPhase->totalMoles())) {
1059  plogf("We have an inconsistency in total moles, %g %g\n",
1060  sumMoles, pubPhase->totalMoles());
1061  exit(EXIT_FAILURE);
1062  }
1063 
1064  }
1065 
1066  pub->m_Iterations = m_VCount->Its;
1068 
1069  return VCS_SUCCESS;
1070 }
1071 /*****************************************************************************/
1072 
1073 // Initialize the internal counters
1074 /*
1075  * Initialize the internal counters containing the subroutine call
1076  * values and times spent in the subroutines.
1077  *
1078  * ifunc = 0 Initialize only those counters appropriate for the top of
1079  * vcs_solve_TP().
1080  * = 1 Initialize all counters.
1081  */
1083 {
1084  m_VCount->Its = 0;
1085  m_VCount->Basis_Opts = 0;
1086  m_VCount->Time_vcs_TP = 0.0;
1087  m_VCount->Time_basopt = 0.0;
1088  if (ifunc) {
1089  m_VCount->T_Its = 0;
1090  m_VCount->T_Basis_Opts = 0;
1091  m_VCount->T_Calls_Inest = 0;
1092  m_VCount->T_Calls_vcs_TP = 0;
1093  m_VCount->T_Time_vcs_TP = 0.0;
1094  m_VCount->T_Time_basopt = 0.0;
1095  m_VCount->T_Time_inest = 0.0;
1096  m_VCount->T_Time_vcs = 0.0;
1097  }
1098 }
1099 /**************************************************************************/
1100 
1101 // Calculation of the total volume and the partial molar volumes
1102 /*
1103  * This function calculates the partial molar volume
1104  * for all species, kspec, in the thermo problem
1105  * at the temperature TKelvin and pressure, Pres, pres is in atm.
1106  * And, it calculates the total volume of the combined system.
1107  *
1108  * Input
1109  * ---------------
1110  * @param tkelvin Temperature in kelvin()
1111  * @param pres Pressure in Pascal
1112  * @param w w[] is thevector containing the current mole numbers
1113  * in units of kmol.
1114  *
1115  * Output
1116  * ----------------
1117  * @param volPM[] For species in all phase, the entries are the
1118  * partial molar volumes units of M**3 / kmol.
1119  *
1120  * @return The return value is the total volume of
1121  * the entire system in units of m**3.
1122  */
1123 double VCS_SOLVE::vcs_VolTotal(const double tkelvin, const double pres,
1124  const double w[], double volPM[])
1125 {
1126  double VolTot = 0.0;
1127  for (size_t iphase = 0; iphase < m_numPhases; iphase++) {
1128  vcs_VolPhase* Vphase = m_VolPhaseList[iphase];
1129  Vphase->setState_TP(tkelvin, pres);
1130  Vphase->setMolesFromVCS(VCS_STATECALC_OLD, w);
1131  double Volp = Vphase->sendToVCS_VolPM(volPM);
1132  VolTot += Volp;
1133  }
1134  return VolTot;
1135 }
1136 /****************************************************************************/
1137 
1138 
1139 }