Cantera  2.0
vcs_TP.cpp
3 #include "vcs_species_thermo.h"
5
6 #include <cstdio>
7 #include <cstdlib>
8 #include <cmath>
9
10 namespace VCSnonideal
11 {
12
13 /*****************************************************************************/
14 /*****************************************************************************/
15 /*****************************************************************************/
16
17 int VCS_SOLVE::vcs_TP(int ipr, int ip1, int maxit, double T_arg, double pres_arg)
18
19 /**************************************************************************
20  *
21  * vcs_TP:
22  *
23  * Solve an equilibrium problem at a particular fixed temperature
24  * and pressure
25  *
26  * ipr = 1 -> Print results to standard output
27  * 0 -> don't report on anything
28  * ip1 = 1 -> Print intermediate results.
29  * maxit -> Maximum number of iterations for the algorithm
30  * T = Temperature (Kelvin)
31  * pres = Pressure (pascal)
32  *
33  * Return Codes
34  * ------------------
35  * 0 = Equilibrium Achieved
36  * 1 = Range space error encountered. The element abundance criteria are
37  * only partially satisfied. Specifically, the first NC= (number of
38  * components) conditions are satisfied. However, the full NE
39  * (number of elements) conditions are not satisfied. The equilibrium
40  * condition is returned.
41  * -1 = Maximum number of iterations is exceeded. Convergence was not
42  * found.
43  ***************************************************************************/
44 {
45  int retn, iconv;
46  /*
47  * Store the temperature and pressure in the private global variables
48  */
49  m_temperature = T_arg;
50  m_pressurePA = pres_arg;
51  /*
52  * Evaluate the standard state free energies
53  * at the current temperatures and pressures.
54  */
55  iconv = vcs_evalSS_TP(ipr, ip1, m_temperature, pres_arg);
56
57  /*
58  * Prepare the problem data:
59  * ->nondimensionalize the free energies using
60  * the divisor, R * T
61  */
62  vcs_nondim_TP();
63  /*
64  * Prep the fe field
65  */
66  vcs_fePrep_TP();
67  /*
68  * Decide whether we need an initial estimate of the solution
69  * If so, go get one. If not, then
70  */
71  if (m_doEstimateEquil) {
72  retn = vcs_inest_TP();
73  if (retn != VCS_SUCCESS) {
74  plogf("vcs_inest_TP returned a failure flag\n");
75  }
76  }
77  /*
78  * Solve the problem at a fixed Temperature and Pressure
79  * (all information concerning Temperature and Pressure has already
80  * been derived. The free energies are now in dimensionless form.)
81  */
82  iconv = vcs_solve_TP(ipr, ip1, maxit);
83
84  /*
85  * Redimensionalize the free energies using
86  * the reverse of vcs_nondim to add back units.
87  */
88  vcs_redim_TP();
89  /*
90  * Return the convergence success flag.
91  */
92  return iconv;
93 }
94 /*****************************************************************************/
95 /*****************************************************************************/
96 /*****************************************************************************/
97 /*ARGSUSED*/
98 int VCS_SOLVE::vcs_evalSS_TP(int ipr, int ip1, double Temp, double pres)
99
100 /**************************************************************************
101  *
102  * vcs_evalSS_TP:
103  *
104  * IPR = 1 -> Print results to standard output
105  * 0 -> don't report on anything
106  * IP1 = 1 -> Print intermediate results.
107  * T = Temperature (Kelvin)
108  * Pres = Pressure (Pascal)
109  *
110  * Evaluate the standard state free energies at the current temperature
111  * and pressure. Ideal gas pressure contribution is added in here.
112  *
113  ***************************************************************************/
114 {
115  // int i;
116  //double R;
117  /*
118  * At this level of the program, we are still using values
119  * for the free energies that have units.
120  */
121  // R = vcsUtil_gasConstant(m_VCS_UnitsFormat);
122
123  /*
124  * We need to special case VCS_UNITS_UNITLESS, here.
125  * cpc_ts_GStar_calc() returns units of Kelvin. Also, the temperature
126  * comes into play in calculating the ideal equation of state
127  * contributions, and other equations of state also. Therefore,
128  * we will emulate the VCS_UNITS_KELVIN case, here by changing
129  * the initial gibbs free energy units to Kelvin before feeding
130  * them to the cpc_ts_GStar_calc() routine. Then, we will revert
131  * them back to unitless at the end of this routine.
132  */
133
134
135  /*
136  * Loop over the species calculating the standard state Gibbs free
137  * energies. -> These are energies that only depend upon the Temperature
138  * and possibly on the pressure (i.e., ideal gas, etc).
139  */
140  // HKM -> We can change this to looks over phases, calling the vcs_VolPhase
141  // object. Working to get rid of VCS_SPECIES_THERMO object
142  //for (i = 0; i < m_numSpeciesTot; ++i) {
143  // VCS_SPECIES_THERMO *spt = SpeciesThermo[i];
144  // ff[i] = R * spt->GStar_R_calc(i, Temp, pres);
145  //}
146
147  for (size_t iph = 0; iph < m_numPhases; iph++) {
148  vcs_VolPhase* vph = m_VolPhaseList[iph];
151  }
152
153  if (m_VCS_UnitsFormat == VCS_UNITS_UNITLESS) {
154  for (size_t i = 0; i < m_numSpeciesTot; ++i) {
155  m_SSfeSpecies[i] /= Temp;
156  }
157  }
158  return VCS_SUCCESS;
159 } /***************************************************************************/
160
161 /*****************************************************************************/
162 /*****************************************************************************/
163 /*****************************************************************************/
164
166
167 /**************************************************************************
168  *
169  *
170  ***************************************************************************/
171 {
172  for (size_t i = 0; i < m_numSpeciesTot; ++i) {
173  /*
174  * For single species phases, initialize the chemical
175  * potential with the value of the standard state chemical
176  * potential. This value doesn't change during the calculation
177  */
178  if (m_SSPhase[i]) {
181  }
182  }
183 } /* vcs_fePrep_TP() ********************************************************/
184
185 }
186