Cantera  2.0
vcs_prep.cpp
Go to the documentation of this file.
1 /**
2  * @file vcs_prep.cpp
3  * This file contains some prepatory functions.
4  */
5 
6 /*
7  * Copyright (2005) Sandia Corporation. Under the terms of
8  * Contract DE-AC04-94AL85000 with Sandia Corporation, the
9  * U.S. Government retains certain rights in this software.
10  */
11 
14 #include "cantera/equil/vcs_prob.h"
16 #include "vcs_SpeciesProperties.h"
17 
18 #include <cstdio>
19 #include <cstdlib>
20 #include <cmath>
21 
22 namespace VCSnonideal
23 {
24 
25 
26 // Calculate the status of single species phases.
28 {
29  size_t iph;
30  vcs_VolPhase* Vphase;
31 
32  std::vector<int> numPhSpecies(m_numPhases, 0);
33 
34  for (size_t kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
35  numPhSpecies[m_phaseID[kspec]]++;
36  }
37  /*
38  * Handle the special case of a single species in a phase that
39  * has been earmarked as a multispecies phase.
40  * Treat that species as a single-species phase
41  */
42  for (iph = 0; iph < m_numPhases; iph++) {
43  Vphase = m_VolPhaseList[iph];
44  Vphase->m_singleSpecies = false;
45  if (TPhInertMoles[iph] > 0.0) {
46  Vphase->setExistence(2);
47  }
48  if (numPhSpecies[iph] <= 1) {
49  if (TPhInertMoles[iph] == 0.0) {
50  Vphase->m_singleSpecies = true;
51  }
52  }
53  }
54 
55  /*
56  * Fill in some useful arrays here that have to do with the
57  * static information concerning the phase ID of species.
58  * SSPhase = Boolean indicating whether a species is in a
59  * single species phase or not.
60  */
61  for (size_t kspec = 0; kspec < m_numSpeciesTot; kspec++) {
62  iph = m_phaseID[kspec];
63  Vphase = m_VolPhaseList[iph];
64  if (Vphase->m_singleSpecies) {
65  m_SSPhase[kspec] = true;
66  } else {
67  m_SSPhase[kspec] = false;
68  }
69  }
70 }
71 /*****************************************************************************/
72 
73 // This routine is mostly concerned with changing the private data
74 // to be consistent with what's needed for solution. It is called one
75 // time for each new problem structure definition.
76 /*
77  * This routine is always followed by vcs_prep(). Therefore, tasks
78  * that need to be done for every call to vcsc() should be placed in
79  * vcs_prep() and not in this routine.
80  *
81  * The problem structure refers to:
82  *
83  * the number and identity of the species.
84  * the formula matrix and thus the number of components.
85  * the number and identity of the phases.
86  * the equation of state
87  * the method and parameters for determining the standard state
88  * The method and parameters for determining the activity coefficients.
89  *
90  * Tasks:
91  * 0) Fill in the SSPhase[] array.
92  * 1) Check to see if any multispecies phases actually have only one
93  * species in that phase. If true, reassign that phase and species
94  * to be a single-species phase.
95  * 2) Determine the number of components in the problem if not already
96  * done so. During this process the order of the species is changed
97  * in the private data structure. All references to the species
98  * properties must employ the ind[] index vector.
99  *
100  * @param printLvl Print level of the routine
101  *
102  * @return the return code
103  * VCS_SUCCESS = everything went OK
104  *
105  */
107 {
108  size_t kspec, i;
109  int retn = VCS_SUCCESS;
110  double pres, test;
111  double* aw, *sa, *sm, *ss;
112  bool modifiedSoln = false;
113  bool conv;
114 
115  m_debug_print_lvl = printLvl;
116 
117 
118  /*
119  * Calculate the Single Species status of phases
120  * Also calculate the number of species per phase
121  */
122  vcs_SSPhase();
123 
124  /*
125  * Set an initial estimate for the number of noncomponent species
126  * equal to nspecies - nelements. This may be changed below
127  */
129  m_numRxnTot = 0;
130  } else {
132  }
135  for (i = 0; i < m_numRxnRdc; ++i) {
137  }
138 
139  for (kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
140  size_t pID = m_phaseID[kspec];
141  size_t spPhIndex = m_speciesLocalPhaseIndex[kspec];
142  vcs_VolPhase* vPhase = m_VolPhaseList[pID];
143  vcs_SpeciesProperties* spProp = vPhase->speciesProperty(spPhIndex);
144  double sz = 0.0;
145  size_t eSize = spProp->FormulaMatrixCol.size();
146  for (size_t e = 0; e < eSize; e++) {
147  sz += fabs(spProp->FormulaMatrixCol[e]);
148  }
149  if (sz > 0.0) {
150  m_spSize[kspec] = sz;
151  } else {
152  m_spSize[kspec] = 1.0;
153  }
154  }
155 
156  /* ***************************************************** */
157  /* **** DETERMINE THE NUMBER OF COMPONENTS ************* */
158  /* ***************************************************** */
159 
160  /*
161  * Obtain a valid estimate of the mole fraction. This will
162  * be used as an initial ordering vector for prioritizing
163  * which species are defined as components.
164  *
165  * If a mole number estimate was supplied from the
166  * input file, use that mole number estimate.
167  *
168  * If a solution estimate wasn't supplied from the input file,
169  * supply an initial estimate for the mole fractions
170  * based on the relative reverse ordering of the
171  * chemical potentials.
172  *
173  * For voltage unknowns, set these to zero for the moment.
174  */
175  test = -1.0e-10;
176  if (m_doEstimateEquil < 0) {
177  double sum = 0.0;
178  for (kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
180  sum += fabs(m_molNumSpecies_old[kspec]);
181  }
182  }
183  if (fabs(sum) < 1.0E-6) {
184  modifiedSoln = true;
185  if (m_pressurePA <= 0.0) {
186  pres = 1.01325E5;
187  } else {
188  pres = m_pressurePA;
189  }
190  retn = vcs_evalSS_TP(0, 0, m_temperature, pres);
191  for (kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
193  m_molNumSpecies_old[kspec] = - m_SSfeSpecies[kspec];
194  } else {
195  m_molNumSpecies_old[kspec] = 0.0;
196  }
197  }
198  }
199  test = -1.0e20;
200  }
201 
202  /*
203  * NC = number of components is in the vcs.h common block
204  * This call to BASOPT doesn't calculate the stoichiometric
205  * reaction matrix.
206  */
207  std::vector<double> awSpace(m_numSpeciesTot + (m_numElemConstraints + 2)*(m_numElemConstraints), 0.0);
208  aw = VCS_DATA_PTR(awSpace);
209  if (aw == NULL) {
210  plogf("vcs_prep_oneTime: failed to get memory: global bailout\n");
211  return VCS_NOMEMORY;
212  }
213  sa = aw + m_numSpeciesTot;
214  sm = sa + m_numElemConstraints;
215  ss = sm + (m_numElemConstraints)*(m_numElemConstraints);
216  retn = vcs_basopt(true, aw, sa, sm, ss, test, &conv);
217  if (retn != VCS_SUCCESS) {
218  plogf("vcs_prep_oneTime:");
219  plogf(" Determination of number of components failed: %d\n",
220  retn);
221  plogf(" Global Bailout!\n");
222  return retn;
223  }
224 
225  if (m_numSpeciesTot >= m_numComponents) {
226  m_numRxnTot = m_numRxnRdc = m_numSpeciesTot - m_numComponents;
227  for (i = 0; i < m_numRxnRdc; ++i) {
228  m_indexRxnToSpecies[i] = m_numComponents + i;
229  }
230  } else {
231  m_numRxnTot = m_numRxnRdc = 0;
232  }
233 
234  /*
235  * The elements might need to be rearranged.
236  */
237  awSpace.resize(m_numElemConstraints + (m_numElemConstraints + 2)*(m_numElemConstraints), 0.0);
238  aw = VCS_DATA_PTR(awSpace);
239  sa = aw + m_numElemConstraints;
240  sm = sa + m_numElemConstraints;
241  ss = sm + (m_numElemConstraints)*(m_numElemConstraints);
242  retn = vcs_elem_rearrange(aw, sa, sm, ss);
243  if (retn != VCS_SUCCESS) {
244  plogf("vcs_prep_oneTime:");
245  plogf(" Determination of element reordering failed: %d\n",
246  retn);
247  plogf(" Global Bailout!\n");
248  return retn;
249  }
250 
251  // If we mucked up the solution unknowns because they were all
252  // zero to start with, set them back to zero here
253  if (modifiedSoln) {
254  for (kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
255  m_molNumSpecies_old[kspec] = 0.0;
256  }
257  }
258  return VCS_SUCCESS;
259 }
260 
261 /*****************************************************************************/
262 
263 // Prepare the object for re-solution
264 /*
265  * This routine is mostly concerned with changing the private data
266  * to be consistent with that needed for solution. It is called for
267  * every invocation of the vcs_solve() except for the cleanup invocation.
268  *
269  * Tasks:
270  * 1) Initialization of arrays to zero.
271  * 2) Calculate total number of moles in all phases
272  *
273  * return code
274  * VCS_SUCCESS = everything went OK
275  * VCS_PUB_BAD = There is an irreconcilable difference in the
276  * public data structure from when the problem was
277  * initially set up.
278  */
280 {
281  /*
282  * Initialize various arrays in the data to zero
283  */
284  vcs_vdzero(m_feSpecies_old, m_numSpeciesTot);
285  vcs_vdzero(m_feSpecies_new, m_numSpeciesTot);
287  vcs_dzero(&(m_deltaMolNumPhase[0][0]), m_numSpeciesTot * m_numPhases);
288  vcs_izero(&(m_phaseParticipation[0][0]), m_numSpeciesTot * m_numPhases);
289  vcs_dzero(VCS_DATA_PTR(m_deltaPhaseMoles), m_numPhases);
290  vcs_dzero(VCS_DATA_PTR(m_tPhaseMoles_new), m_numPhases);
291  /*
292  * Calculate the total number of moles in all phases.
293  */
294  vcs_tmoles();
295  return VCS_SUCCESS;
296 }
297 
298 /*****************************************************************************/
299 
300 // In this routine, we check for things that will cause the algorithm
301 // to fail.
302 /*
303  * We check to see if the problem is well posed. If it is not, we return
304  * false and print out error conditions.
305  *
306  * Current there is one condition. If all the element abundances are
307  * zero, the algorithm will fail
308  *
309  * @param vprob VCS_PROB pointer to the definition of the equilibrium
310  * problem
311  *
312  * @return If true, the problem is well-posed. If false, the problem
313  * is not well posed.
314  */
316 {
317  double sum = 0.0;
318  for (size_t e = 0; e < vprob->ne; e++) {
319  sum = sum + vprob->gai[e];
320  }
321  if (sum < 1.0E-20) {
322  plogf("vcs_wellPosed: Element abundance is close to zero\n");
323  return false;
324  }
325  return true;
326 }
327 
328 /*****************************************************************************/
329 }