1 /** 2 * @file vcs_solve.h Header file for the internal object that holds the vcs 3 * equilibrium problem (see Class \link Cantera::VCS_SOLVE 4 * VCS_SOLVE\endlink and \ref equilfunctions ). 5 */ 6 7 // This file is part of Cantera. See License.txt in the top-level directory or 8 // at https://cantera.org/license.txt for license and copyright information. 9 10 #ifndef _VCS_SOLVE_H 11 #define _VCS_SOLVE_H 12 13 /* 14 * Index of Symbols 15 * ------------------- 16 * irxn -> refers to the species or rxn between the species and 17 * the components in the problem 18 * k -> refers to the species 19 * j -> refers to the element or component 20 * 21 * ### -> to be eliminated 22 */ 23 24 #include "cantera/base/ct_defs.h" 25 #include "cantera/equil/vcs_defs.h" 26 #include "cantera/equil/vcs_internal.h" 27 #include "cantera/base/Array.h" 28 29 namespace Cantera 30 { 31 32 class vcs_VolPhase; 33 class VCS_SPECIES_THERMO; 34 class VCS_COUNTERS; 35 class MultiPhase; 36 37 //! This is the main structure used to hold the internal data 38 //! used in vcs_solve_TP(), and to solve TP systems. 39 /*! 40 * The indices of information in this structure may change when the species 41 * basis changes or when phases pop in and out of existence. Both of these 42 * operations change the species ordering. 43 */ 44 class VCS_SOLVE 45 { 46 public: 47 //! Initialize the sizes within the VCS_SOLVE object 48 /*! 49 * This resizes all of the internal arrays within the object. 50 */ 51 VCS_SOLVE(MultiPhase* mphase, int printLvl=0); 52 53 ~VCS_SOLVE(); 54 55 //! Solve an equilibrium problem 56 /*! 57 * This is the main interface routine to the equilibrium solver 58 * 59 * @param ipr Printing of results 60 * ipr = 1 -> Print problem statement and final results to 61 * standard output 62 * 0 -> don't report on anything 63 * @param ip1 Printing of intermediate results 64 * IP1 = 1 -> Print intermediate results. 65 * @param maxit Maximum number of iterations for the algorithm 66 * @return nonzero value: failure to solve the problem at hand. zero : 67 * success 68 */ 69 int vcs(int ipr, int ip1, int maxit); 70 71 //! Main routine that solves for equilibrium at constant T and P using a 72 //! variant of the VCS method 73 /*! 74 * This is the main routine that solves for equilibrium at constant T and P 75 * using a variant of the VCS method. Nonideal phases can be accommodated as 76 * well. 77 * 78 * Any number of single-species phases and multi-species phases can be 79 * handled by the present version. 80 * 81 * @param print_lvl 1 -> Print results to standard output; 82 * 0 -> don't report on anything 83 * @param printDetails 1 -> Print intermediate results. 84 * @param maxit Maximum number of iterations for the algorithm 85 * @return 86 * * 0 = Equilibrium Achieved 87 * * 1 = Range space error encountered. The element abundance criteria 88 * are only partially satisfied. Specifically, the first NC= (number 89 * of components) conditions are satisfied. However, the full NE 90 * (number of elements) conditions are not satisfied. The 91 * equilibrium condition is returned. 92 * * -1 = Maximum number of iterations is exceeded. Convergence was 93 * not found. 94 */ 95 int vcs_solve_TP(int print_lvl, int printDetails, int maxit); 96 97 /*! 98 * We make decisions on the initial mole number, and major-minor status 99 * here. We also fix up the total moles in a phase. 100 * 101 * irxn = id of the noncomponent species formation reaction for the 102 * species to be added in. 103 * 104 * The algorithm proceeds to implement these decisions in the previous 105 * position of the species. Then, vcs_switch_pos is called to move the 106 * species into the last active species slot, incrementing the number of 107 * active species at the same time. 108 * 109 * This routine is responsible for the global data manipulation only. 110 */ 111 void vcs_reinsert_deleted(size_t kspec); 112 113 //! Choose the optimum species basis for the calculations 114 /*! 115 * This is done by choosing the species with the largest mole fraction not 116 * currently a linear combination of the previous components. Then, 117 * calculate the stoichiometric coefficient matrix for that basis. 118 * 119 * Rearranges the solution data to put the component data at the 120 * front of the species list. 121 * 122 * Then, calculates `m_stoichCoeffRxnMatrix(jcomp,irxn)` the formation 123 * reactions for all noncomponent species in the mechanism. Also 124 * calculates DNG(I) and DNL(I), the net mole change for each formation 125 * reaction. Also, initializes IR(I) to the default state. 126 * 127 * @param[in] doJustComponents If true, the #m_stoichCoeffRxnMatrix and 128 * #m_deltaMolNumPhase are not calculated. 129 * @param[in] aw Vector of mole fractions which will be used to 130 * construct an optimal basis from. 131 * @param[in] sa Gram-Schmidt orthog work space (nc in length) sa[j] 132 * @param[in] ss Gram-Schmidt orthog work space (nc in length) ss[j] 133 * @param[in] sm QR matrix work space (nc*ne in length) sm[i+j*ne] 134 * @param[in] test This is a small negative number dependent upon whether 135 * an estimate is supplied or not. 136 * @param[out] usedZeroedSpecies If true, then a species with a zero 137 * concentration was used as a component. 138 * The problem may be converged. Or, the 139 * problem may have a range space error and 140 * may not have a proper solution. 141 * @returns VCS_SUCCESS if everything went ok. Returns 142 * VCS_FAILED_CONVERGENCE if there is a problem. 143 * 144 * ### Internal Variables calculated by this routine: 145 * 146 * - #m_numComponents: Number of component species. This routine 147 * calculates the #m_numComponents species. It switches their positions 148 * in the species vector so that they occupy the first #m_numComponents 149 * spots in the species vector. 150 * - `m_stoichCoeffRxnMatrix(jcomp,irxn)` Stoichiometric coefficient 151 * matrix for the reaction mechanism expressed in Reduced Canonical 152 * Form. jcomp refers to the component number, and irxn refers to the 153 * irxn_th non-component species. 154 * - `m_deltaMolNumPhase(iphase,irxn)`: Change in the number of moles in 155 * phase, iphase, due to the noncomponent formation reaction, irxn. 156 * - `m_phaseParticipation(iphase,irxn)`: This is 1 if the phase, iphase, 157 * participates in the formation reaction, irxn, and zero otherwise. 158 */ 159 int vcs_basopt(const bool doJustComponents, double aw[], double sa[], double sm[], 160 double ss[], double test, bool* const usedZeroedSpecies); 161 162 //! Choose a species to test for the next component 163 /*! 164 * We make the choice based on testing (molNum[i] * spSize[i]) for its 165 * maximum value. Preference for single species phases is also made. 166 * 167 * @param molNum Mole number vector 168 * @param j index into molNum[] that indicates where the search will 169 * start from Previous successful components are swapped into the front 170 * of molNum[]. 171 * @param n Length of molNum[] 172 */ 173 size_t vcs_basisOptMax(const double* const molNum, const size_t j, const size_t n); 174 175 //! Evaluate the species category for the indicated species 176 /*! 177 * All evaluations are done using the "old" version of the solution. 178 * 179 * @param kspec Species to be evaluated 180 * @returns the calculated species type 181 */ 182 int vcs_species_type(const size_t kspec) const; 183 184 //! This routine evaluates the species type for all species 185 /*! 186 * - #VCS_SPECIES_MAJOR: Major species 187 * - #VCS_SPECIES_MINOR: Minor species 188 * - #VCS_SPECIES_SMALLMS: The species lies in a multicomponent phase 189 * that exists. Its concentration is currently very low, necessitating 190 * a different method of calculation. 191 * - #VCS_SPECIES_ZEROEDMS: The species lies in a multicomponent phase 192 * which currently doesn't exist. Its concentration is currently zero. 193 * - #VCS_SPECIES_ZEROEDSS: Species lies in a single-species phase which 194 * is currently zeroed out. 195 * - #VCS_SPECIES_DELETED: Species has such a small mole fraction it is 196 * deleted even though its phase may possibly exist. The species is 197 * believed to have such a small mole fraction that it best to throw 198 * the calculation of it out. It will be added back in at the end of 199 * the calculation. 200 * - #VCS_SPECIES_INTERFACIALVOLTAGE: Species refers to an electron in 201 * the metal The unknown is equal to the interfacial voltage drop 202 * across the interface on the SHE (standard hydrogen electrode) scale 203 * (volts). 204 * - #VCS_SPECIES_ZEROEDPHASE: Species lies in a multicomponent phase 205 * that is zeroed atm and will stay deleted due to a choice from a 206 * higher level. These species will formally always have zero mole 207 * numbers in the solution vector. 208 * - #VCS_SPECIES_ACTIVEBUTZERO: The species lies in a multicomponent 209 * phase which currently does exist. Its concentration is currently 210 * identically zero, though the phase exists. Note, this is a temporary 211 * condition that exists at the start of an equilibrium problem. The 212 * species is soon "birthed" or "deleted". 213 * - #VCS_SPECIES_STOICHZERO: The species lies in a multicomponent phase 214 * which currently does exist. Its concentration is currently 215 * identically zero, though the phase exists. This is a permanent 216 * condition due to stoich constraints 217 */ 218 bool vcs_evaluate_speciesType(); 219 220 //! Calculate the dimensionless chemical potentials of all species or 221 //! of certain groups of species, at a fixed temperature and pressure. 222 /*! 223 * We calculate the dimensionless chemical potentials of all species or 224 * certain groups of species here, at a fixed temperature and pressure, for 225 * the input mole vector z[] in the parameter list. Nondimensionalization is 226 * achieved by division by RT. 227 * 228 * Note, for multispecies phases which are currently zeroed out, the 229 * chemical potential is filled out with the standard chemical potential. 230 * 231 * For species in multispecies phases whose concentration is zero, we need 232 * to set the mole fraction to a very low value. Its chemical potential is 233 * then calculated using the VCS_DELETE_MINORSPECIES_CUTOFF concentration 234 * to keep numbers positive. 235 * 236 * For formulas, see vcs_chemPotPhase(). 237 * 238 * Handling of Small Species: 239 * ------------------------------ 240 * As per the discussion above, for small species where the mole fraction 241 * 242 * z(i) < VCS_DELETE_MINORSPECIES_CUTOFF 243 * 244 * The chemical potential is calculated as: 245 * 246 * m_feSpecies(I)(I) = m_SSfeSpecies(I) + ln(ActCoeff[i](VCS_DELETE_MINORSPECIES_CUTOFF)) 247 * 248 * Species in the following categories are treated as "small species" 249 * - #VCS_SPECIES_DELETED 250 * - #VCS_SPECIES_ACTIVEBUTZERO 251 * 252 * For species in multispecies phases which are currently not active, the 253 * treatment is different. These species are in the following species 254 * categories: 255 * - #VCS_SPECIES_ZEROEDMS 256 * - #VCS_SPECIES_ZEROEDPHASE 257 * 258 * For these species, the `ln( ActCoeff[I] X[I])` term is dropped 259 * altogether. The following equation is used: 260 * 261 * m_feSpecies(I) = m_SSfeSpecies(I) 262 * + Charge[I] * Faraday_dim * phasePhi[iphase]; 263 * 264 * Handling of "Species" Representing Interfacial Voltages 265 * --------------------------------------------------------- 266 * 267 * These species have species types of 268 * #VCS_SPECIES_TYPE_INTERFACIALVOLTAGE The chemical potentials for these 269 * "species" refer to electrons in metal electrodes. They have the 270 * following formula 271 * 272 * m_feSpecies(I) = m_SSfeSpecies(I) - F z[I] / RT 273 * 274 * - `F` is Faraday's constant. 275 * - `R` = gas constant 276 * - `T` = temperature 277 * - `V` = potential of the interface = phi_electrode - phi_solution 278 * 279 * For these species, the solution vector unknown, z[I], is V, the phase 280 * voltage, in volts. 281 * 282 * @param ll Determine which group of species gets updated 283 * - `ll = 0`: Calculate for all species 284 * - `ll < 0`: calculate for components and for major non-components 285 * - `ll = 1`: calculate for components and for minor non-components 286 * 287 * @param lbot Restricts the calculation of the chemical potential 288 * to the species between LBOT <= i < LTOP. Usually 289 * LBOT and LTOP will be equal to 0 and MR, respectively. 290 * @param ltop Top value of the loops 291 * 292 * @param stateCalc Determines whether z is old or new or tentative: 293 * - 1: Use the tentative values for the total number of 294 * moles in the phases, i.e., use TG1 instead of TG etc. 295 * - 0: Use the base values of the total number of 296 * moles in each system. 297 * 298 * Also needed: 299 * ff : standard state chemical potentials. These are the 300 * chemical potentials of the standard states at 301 * the same T and P as the solution. 302 * tg : Total Number of moles in the phase. 303 */ 304 void vcs_dfe(const int stateCalc, const int ll, const size_t lbot, const size_t ltop); 305 306 //! This routine uploads the state of the system into all of the 307 //! vcs_VolumePhase objects in the current problem. 308 /*! 309 * @param stateCalc Determines where to get the mole numbers from. 310 * - VCS_STATECALC_OLD -> from m_molNumSpecies_old 311 * - VCS_STATECALC_NEW -> from m_molNumSpecies_new 312 */ 313 void vcs_updateVP(const int stateCalc); 314 315 //! Utility function that evaluates whether a phase can be popped 316 //! into existence 317 /*! 318 * A phase can be popped iff the stoichiometric coefficients for the 319 * component species, whose concentrations will be lowered during the 320 * process, are positive by at least a small degree. 321 * 322 * If one of the phase species is a zeroed component, then the phase can 323 * be popped if the component increases in mole number as the phase moles 324 * are increased. 325 * 326 * @param iphasePop id of the phase, which is currently zeroed, 327 * @returns true if the phase can come into existence and false otherwise. 328 */ 329 bool vcs_popPhasePossible(const size_t iphasePop) const; 330 331 //! Decision as to whether a phase pops back into existence 332 /*! 333 * @param phasePopPhaseIDs Vector containing the phase ids of the phases 334 * that will be popped this step. 335 * @returns the phase id of the phase that pops back into existence. Returns 336 * -1 if there are no phases 337 */ 338 size_t vcs_popPhaseID(std::vector<size_t> &phasePopPhaseIDs); 339 340 //! Calculates the deltas of the reactions due to phases popping 341 //! into existence 342 /*! 343 * Updates #m_deltaMolNumSpecies[irxn] : reaction adjustments, where irxn 344 * refers to the irxn'th species formation reaction. This adjustment is 345 * for species irxn + M, where M is the number of components. 346 * 347 * @param iphasePop Phase id of the phase that will come into existence 348 * @returns an int representing the status of the step 349 * - 0 : normal return 350 * - 1 : A single species phase species has been zeroed out 351 * in this routine. The species is a noncomponent 352 * - 2 : Same as one but, the zeroed species is a component. 353 * - 3 : Nothing was done because the phase couldn't be birthed 354 * because a needed component is zero. 355 */ 356 int vcs_popPhaseRxnStepSizes(const size_t iphasePop); 357 358 //! Calculates formation reaction step sizes. 359 /*! 360 * This is equation 6.4-16, p. 143 in Smith and Missen. 361 * 362 * Output 363 * ------- 364 * m_deltaMolNumSpecies(irxn) : reaction adjustments, where irxn refers to 365 * the irxn'th species formation reaction. 366 * This adjustment is for species irxn + M, 367 * where M is the number of components. 368 * 369 * Special branching occurs sometimes. This causes the component basis 370 * to be reevaluated 371 * 372 * @param forceComponentCalc integer flagging whether a component 373 * recalculation needs to be carried out. 374 * @param kSpecial species number of phase being zeroed. 375 * @returns an int representing which phase may need to be zeroed 376 */ 377 size_t vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial); 378 379 //! Calculates the total number of moles of species in all phases. 380 /*! 381 * Also updates the variable m_totalMolNum and Reconciles Phase existence 382 * flags with total moles in each phase. 383 */ 384 double vcs_tmoles(); 385 386 void check_tmoles() const; 387 388 //! This subroutine calculates reaction free energy changes for 389 //! all noncomponent formation reactions. 390 /*! 391 * Formation reactions are 392 * reactions which create each noncomponent species from the component 393 * species. `m_stoichCoeffRxnMatrix(jcomp,irxn)` are the stoichiometric 394 * coefficients for these reactions. A stoichiometric coefficient of 395 * one is assumed for species irxn in this reaction. 396 * 397 * @param L 398 * - `L < 0`: Calculate reactions corresponding to major noncomponent 399 * and zeroed species only 400 * - `L = 0`: Do all noncomponent reactions, i, between 401 * 0 <= i < irxnl 402 * - `L > 0`: Calculate reactions corresponding to minor noncomponent 403 * and zeroed species only 404 * 405 * @param doDeleted Do deleted species 406 * @param vcsState Calculate deltaG corresponding to either old or new 407 * free energies 408 * @param alterZeroedPhases boolean indicating whether we should 409 * add in a special section for zeroed phases. 410 * 411 * Note we special case one important issue. If the component has zero 412 * moles, then we do not allow deltaG < 0.0 for formation reactions which 413 * would lead to the loss of more of that same component. This dG < 0.0 414 * condition feeds back into the algorithm in several places, and leads 415 * to a infinite loop in at least one case. 416 */ 417 void vcs_deltag(const int L, const bool doDeleted, const int vcsState, 418 const bool alterZeroedPhases = true); 419 420 void vcs_printDeltaG(const int stateCalc); 421 422 //! Swaps the indices for all of the global data for two species, k1 423 //! and k2. 424 /*! 425 * @param ifunc If true, switch the species data and the noncomponent 426 * reaction data. This must be called for a non-component species only. 427 * If false, switch the species data only. Typically, we use this option 428 * when determining the component species and at the end of the 429 * calculation, when we want to return unscrambled results. All rxn data 430 * will be out-of-date. 431 * @param k1 First species index 432 * @param k2 Second species index 433 */ 434 void vcs_switch_pos(const bool ifunc, const size_t k1, const size_t k2); 435 436 //! Main program to test whether a deleted phase should be brought 437 //! back into existence 438 /*! 439 * @param iph Phase id of the deleted phase 440 */ 441 double vcs_phaseStabilityTest(const size_t iph); 442 443 //! Solve an equilibrium problem at a particular fixed temperature 444 //! and pressure 445 /*! 446 * The actual problem statement is assumed to be in the structure already. 447 * This is a wrapper around the solve_TP() function. In this wrapper, we 448 * calculate the standard state Gibbs free energies of the species 449 * and we decide whether to we need to use the initial guess algorithm. 450 * 451 * @param ipr = 1 -> Print results to standard output; 452 * 0 -> don't report on anything 453 * @param ip1 = 1 -> Print intermediate results; 454 * 0 -> Dont print any intermediate results 455 * @param maxit Maximum number of iterations for the algorithm 456 * @param T Value of the Temperature (Kelvin) 457 * @param pres Value of the Pressure 458 * @returns an integer representing the success of the algorithm 459 * * 0 = Equilibrium Achieved 460 * * 1 = Range space error encountered. The element abundance criteria are 461 * only partially satisfied. Specifically, the first NC= (number of 462 * components) conditions are satisfied. However, the full NE (number of 463 * elements) conditions are not satisfied. The equilibrium condition is 464 * returned. 465 * * -1 = Maximum number of iterations is exceeded. Convergence was not 466 * found. 467 */ 468 int vcs_TP(int ipr, int ip1, int maxit, double T, double pres); 469 470 /*! 471 * Evaluate the standard state free energies at the current temperature 472 * and pressure. Ideal gas pressure contribution is added in here. 473 * 474 * @param ipr 1 -> Print results to standard output; 0 -> don't report 475 * on anything 476 * @param ip1 1 -> Print intermediate results; 0 -> don't. 477 * @param Temp Temperature (Kelvin) 478 * @param pres Pressure (Pascal) 479 */ 480 int vcs_evalSS_TP(int ipr, int ip1, double Temp, double pres); 481 482 //! Initialize the chemical potential of single species phases 483 /*! 484 * For single species phases, initialize the chemical potential with the 485 * value of the standard state chemical potential. This value doesn't 486 * change during the calculation 487 */ 488 void vcs_fePrep_TP(); 489 490 //! Calculation of the total volume and the partial molar volumes 491 /*! 492 * This function calculates the partial molar volume for all species, kspec, 493 * in the thermo problem at the temperature TKelvin and pressure, Pres, pres 494 * is in atm. And, it calculates the total volume of the combined system. 495 * 496 * @param[in] tkelvin Temperature in kelvin() 497 * @param[in] pres Pressure in Pascal 498 * @param[in] w w[] is the vector containing the current mole 499 * numbers in units of kmol. 500 * @param[out] volPM[] For species in all phase, the entries are the 501 * partial molar volumes units of M**3 / kmol. 502 * @returns the total volume of the entire system in units of m**3. 503 */ 504 double vcs_VolTotal(const double tkelvin, const double pres, 505 const double w[], double volPM[]); 506 507 //! This routine is mostly concerned with changing the private data to be 508 //! consistent with what's needed for solution. It is called one time for 509 //! each new problem structure definition. 510 /*! 511 * The problem structure refers to: 512 * 513 * - the number and identity of the species. 514 * - the formula matrix and thus the number of components. 515 * - the number and identity of the phases. 516 * - the equation of state 517 * - the method and parameters for determining the standard state 518 * - The method and parameters for determining the activity coefficients. 519 * 520 * Tasks: 521 * 1. Fill in the SSPhase[] array. 522 * 2. Check to see if any multispecies phases actually have only one 523 * species in that phase. If true, reassign that phase and species 524 * to be a single-species phase. 525 * 3. Determine the number of components in the problem if not already 526 * done so. During this process the order of the species is changed 527 * in the private data structure. All references to the species 528 * properties must employ the ind[] index vector. 529 * 4. Initialization of arrays to zero. 530 * 5. Check to see if the problem is well posed (If all the element 531 * abundances are zero, the algorithm will fail) 532 * 533 * @param printLvl Print level of the routine 534 * @return 535 * VCS_SUCCESS = everything went OK; 536 * VCS_PUB_BAD = There is an irreconcilable difference in the 537 * public data structure from when the problem was 538 * initially set up. 539 */ 540 int vcs_prep(int printLvl); 541 542 //! Rearrange the constraint equations represented by the Formula 543 //! Matrix so that the operational ones are in the front 544 /*! 545 * This subroutine handles the rearrangement of the constraint equations 546 * represented by the Formula Matrix. Rearrangement is only necessary when 547 * the number of components is less than the number of elements. For this 548 * case, some constraints can never be satisfied exactly, because the 549 * range space represented by the Formula Matrix of the components can't 550 * span the extra space. These constraints, which are out of the range 551 * space of the component Formula matrix entries, are migrated to the back 552 * of the Formula matrix. 553 * 554 * A prototypical example is an extra element column in FormulaMatrix[], 555 * which is identically zero. For example, let's say that argon is has an 556 * element column in FormulaMatrix[], but no species in the mechanism 557 * actually contains argon. Then, nc < ne. Also, without perturbation of 558 * FormulaMatrix[] vcs_basopt[] would produce a zero pivot because the 559 * matrix would be singular (unless the argon element column was already 560 * the last column of FormulaMatrix[]. 561 * 562 * This routine borrows heavily from vcs_basopt's algorithm. It finds nc 563 * constraints which span the range space of the Component Formula matrix, 564 * and assigns them as the first nc components in the formula matrix. This 565 * guarantees that vcs_basopt[] has a nonsingular matrix to invert. 566 * 567 * Other Variables 568 * @param aw Mole fraction work space (ne in length) 569 * @param sa Gram-Schmidt orthog work space (ne in length) 570 * @param sm QR matrix work space (ne*ne in length) 571 * @param ss Gram-Schmidt orthog work space (ne in length) 572 */ 573 int vcs_elem_rearrange(double* const aw, double* const sa, 574 double* const sm, double* const ss); 575 576 //! Swaps the indices for all of the global data for two elements, ipos 577 //! and jpos. 578 /*! 579 * This function knows all of the element information with VCS_SOLVE, and 580 * can therefore switch element positions 581 * 582 * @param ipos first global element index 583 * @param jpos second global element index 584 */ 585 void vcs_switch_elem_pos(size_t ipos, size_t jpos); 586 587 //! Calculates the diagonal contribution to the Hessian due to 588 //! the dependence of the activity coefficients on the mole numbers. 589 /*! 590 * (See framemaker notes, Eqn. 20 - VCS Equations document) 591 * 592 * We allow the diagonal to be increased positively to any degree. 593 * We allow the diagonal to be decreased to 1/3 of the ideal solution 594 * value, but no more -> it must remain positive. 595 */ 596 double vcs_Hessian_diag_adj(size_t irxn, double hessianDiag_Ideal); 597 598 //! Calculates the diagonal contribution to the Hessian due to 599 //! the dependence of the activity coefficients on the mole numbers. 600 /*! 601 * (See framemaker notes, Eqn. 20 - VCS Equations document) 602 */ 603 double vcs_Hessian_actCoeff_diag(size_t irxn); 604 605 //! Recalculate all of the activity coefficients in all of the phases 606 //! based on input mole numbers 607 /* 608 * @param moleSpeciesVCS kmol of species to be used in the update. 609 * 610 * NOTE: This routine needs to be regulated. 611 */ 612 void vcs_CalcLnActCoeffJac(const double* const moleSpeciesVCS); 613 614 //! Print out a report on the state of the equilibrium problem to 615 //! standard output. 616 /*! 617 * @param iconv Indicator of convergence, to be printed out in the report: 618 * - 0 converged 619 * - 1 range space error 620 * - -1 not converged 621 */ 622 int vcs_report(int iconv); 623 624 //! Computes the current elemental abundances vector 625 /*! 626 * Computes the elemental abundances vector, m_elemAbundances[], and stores 627 * it back into the global structure 628 */ 629 void vcs_elab(); 630 631 /*! 632 * Checks to see if the element abundances are in compliance. If they are, 633 * then TRUE is returned. If not, FALSE is returned. Note the number of 634 * constraints checked is usually equal to the number of components in the 635 * problem. This routine can check satisfaction of all of the constraints 636 * in the problem, which is equal to ne. However, the solver can't fix 637 * breakage of constraints above nc, because that nc is the range space by 638 * definition. Satisfaction of extra constraints would have had to occur 639 * in the problem specification. 640 * 641 * The constraints should be broken up into 2 sections. If a constraint 642 * involves a formula matrix with positive and negative signs, and eaSet = 643 * 0.0, then you can't expect that the sum will be zero. There may be 644 * roundoff that inhibits this. However, if the formula matrix is all of 645 * one sign, then this requires that all species with nonzero entries in 646 * the formula matrix be identically zero. We put this into the logic 647 * below. 648 * 649 * @param ibound 1: Checks constraints up to the number of elements; 650 * 0: Checks constraints up to the number of components. 651 */ 652 bool vcs_elabcheck(int ibound); 653 654 /*! 655 * Computes the elemental abundances vector for a single phase, 656 * elemAbundPhase[], and returns it through the argument list. The mole 657 * numbers of species are taken from the current value in 658 * m_molNumSpecies_old[]. 659 */ 660 void vcs_elabPhase(size_t iphase, double* const elemAbundPhase); 661 662 /*! 663 * This subroutine corrects for element abundances. At the end of the 664 * subroutine, the total moles in all phases are recalculated again, 665 * because we have changed the number of moles in this routine. 666 * 667 * @param aa temporary work vector, length ne*ne 668 * @param x temporary work vector, length ne 669 * 670 * @returns 671 * - 0 = Nothing of significance happened, Element abundances were and 672 * still are good. 673 * - 1 = The solution changed significantly; The element abundances are 674 * now good. 675 * - 2 = The solution changed significantly, The element abundances are 676 * still bad. 677 * - 3 = The solution changed significantly, The element abundances are 678 * still bad and a component species got zeroed out. 679 * 680 * Internal data to be worked on:: 681 * - ga Current element abundances 682 * - m_elemAbundancesGoal Required elemental abundances 683 * - m_molNumSpecies_old Current mole number of species. 684 * - m_formulaMatrix[][] Formula matrix of the species 685 * - ne Number of elements 686 * - nc Number of components. 687 * 688 * NOTES: 689 * This routine is turning out to be very problematic. There are 690 * lots of special cases and problems with zeroing out species. 691 * 692 * Still need to check out when we do loops over nc vs. ne. 693 */ 694 int vcs_elcorr(double aa[], double x[]); 695 696 //! Create an initial estimate of the solution to the thermodynamic 697 //! equilibrium problem. 698 /*! 699 * @return value indicates success: 700 * - 0: successful initial guess 701 * - -1: Unsuccessful initial guess; the elemental abundances aren't 702 * satisfied. 703 */ 704 int vcs_inest_TP(); 705 706 //! Estimate the initial mole numbers by constrained linear programming 707 /*! 708 * This is done by running each reaction as far forward or backward as 709 * possible, subject to the constraint that all mole numbers remain non- 710 * negative. Reactions for which \f$ \Delta \mu^0 \f$ are positive are run 711 * in reverse, and ones for which it is negative are run in the forward 712 * direction. The end result is equivalent to solving the linear 713 * programming problem of minimizing the linear Gibbs function subject to 714 * the element and non-negativity constraints. 715 */ 716 int vcs_setMolesLinProg(); 717 718 //! Calculate the total dimensionless Gibbs free energy 719 /*! 720 * Inert species are handled as if they had a standard free energy of 721 * zero. Note, for this algorithm this function should be monotonically 722 * decreasing. 723 */ 724 double vcs_Total_Gibbs(double* w, double* fe, double* tPhMoles); 725 726 //! Calculate the total dimensionless Gibbs free energy of a single phase 727 /*! 728 * Inert species are handled as if they had a standard free energy of 729 * zero and if they obeyed ideal solution/gas theory. 730 * 731 * @param iphase ID of the phase 732 * @param w Species mole number vector for all species 733 * @param fe vector of partial molar free energies of all of the 734 * species 735 */ 736 double vcs_GibbsPhase(size_t iphase, const double* const w, 737 const double* const fe); 738 739 //! Transfer the results of the equilibrium calculation back from VCS_SOLVE 740 void vcs_prob_update(); 741 742 //! Fully specify the problem to be solved 743 void vcs_prob_specifyFully(); 744 745 private: 746 //! Zero out the concentration of a species. 747 /*! 748 * Make sure to conserve elements and keep track of the total moles in all 749 * phases. 750 * - w[] 751 * - m_tPhaseMoles_old[] 752 * 753 * @param kspec Species index 754 * @return: 755 * 1: succeeded 756 * 0: failed. 757 */ 758 int vcs_zero_species(const size_t kspec); 759 760 //! Change a single species from active to inactive status 761 /*! 762 * Rearrange data when species is added or removed. The Lth species is 763 * moved to the back of the species vector. The back of the species 764 * vector is indicated by the value of MR, the current number of 765 * active species in the mechanism. 766 * 767 * @param kspec Species Index 768 * @return 769 * Returns 0 unless. 770 * The return is 1 when the current number of 771 * noncomponent species is equal to zero. A recheck of deleted species 772 * is carried out in the main code. 773 */ 774 int vcs_delete_species(const size_t kspec); 775 776 //! This routine handles the bookkeeping involved with the deletion of 777 //! multiphase phases from the problem. 778 /*! 779 * When they are deleted, all of their species become active species, even 780 * though their mole numbers are set to zero. The routine does not make the 781 * decision to eliminate multiphases. 782 * 783 * Note, species in phases with zero mole numbers are still considered 784 * active. Whether the phase pops back into existence or not is checked as 785 * part of the main iteration loop. 786 * 787 * @param iph Phase to be deleted 788 * @returns whether the operation was successful or not 789 */ 790 bool vcs_delete_multiphase(const size_t iph); 791 792 //! Change the concentration of a species by delta moles. 793 /*! 794 * Make sure to conserve elements and keep track of the total kmoles in 795 * all phases. 796 * 797 * @param kspec The species index 798 * @param delta_ptr pointer to the delta for the species. This may 799 * change during the calculation 800 * @return 801 * 1: succeeded without change of dx 802 * 0: Had to adjust dx, perhaps to zero, in order to do the delta. 803 */ 804 int delta_species(const size_t kspec, double* const delta_ptr); 805 806 //! Provide an estimate for the deleted species in phases that are not 807 //! zeroed out 808 /*! 809 * Try to add back in all deleted species. An estimate of the kmol numbers 810 * are obtained and the species is added back into the equation system, into 811 * the old state vector. 812 * 813 * This routine is called at the end of the calculation, just before 814 * returning to the user. 815 */ 816 size_t vcs_add_all_deleted(); 817 818 //! Recheck deleted species in multispecies phases. 819 /*! 820 * We are checking the equation: 821 * 822 * sum_u = sum_j_comp [ sigma_i_j * u_j ] 823 * = u_i_O + log((AC_i * W_i)/m_tPhaseMoles_old) 824 * 825 * by first evaluating: 826 * 827 * DG_i_O = u_i_O - sum_u. 828 * 829 * Then, if TL is zero, the phase pops into existence if DG_i_O < 0. Also, 830 * if the phase exists, then we check to see if the species can have a mole 831 * number larger than VCS_DELETE_SPECIES_CUTOFF (default value = 1.0E-32). 832 */ 833 int vcs_recheck_deleted(); 834 835 //! Minor species alternative calculation 836 /*! 837 * This is based upon the following approximation: The mole fraction 838 * changes due to these reactions don't affect the mole numbers of the 839 * component species. Therefore the following approximation is valid for 840 * a small component of an ideal phase: 841 * 842 * 0 = m_deltaGRxn_old(I) + log(molNum_new(I)/molNum_old(I)) 843 * 844 * `m_deltaGRxn_old` contains the contribution from 845 * 846 * m_feSpecies_old(I) = 847 * m_SSfeSpecies(I) + 848 * log(ActCoeff[i] * molNum_old(I) / m_tPhaseMoles_old(iph)) 849 * Thus, 850 * 851 * molNum_new(I)= molNum_old(I) * EXP(-m_deltaGRxn_old(I)) 852 * 853 * Most of this section is mainly restricting the update to reasonable 854 * values. We restrict the update a factor of 1.0E10 up and 1.0E-10 down 855 * because we run into trouble with the addition operator due to roundoff if 856 * we go larger than ~1.0E15. Roundoff will then sometimes produce zero mole 857 * fractions. 858 * 859 * Note: This routine was generalized to incorporate nonideal phases and 860 * phases on the molality basis 861 * 862 * @param[in] kspec The current species and corresponding formation 863 * reaction number. 864 * @param[in] irxn The current species and corresponding formation 865 * reaction number. 866 * @param[out] do_delete: BOOLEAN which if true on return, then we 867 * branch to the section that deletes a species 868 * from the current set of active species. 869 */ 870 double vcs_minor_alt_calc(size_t kspec, size_t irxn, bool* do_delete, 871 char* ANOTE=0) const; 872 873 //! This routine optimizes the minimization of the total Gibbs free energy 874 //! by making sure the slope of the following functional stays negative: 875 /*! 876 * The slope of the following functional is equivalent to the slope of the 877 * total Gibbs free energy of the system: 878 * 879 * d_Gibbs/ds = sum_k( m_deltaGRxn * m_deltaMolNumSpecies[k] ) 880 * 881 * along the current direction m_deltaMolNumSpecies[], by choosing a value, al: (0<al<1) 882 * such that the a parabola approximation to Gibbs(al) fit to the 883 * end points al = 0 and al = 1 is minimized. 884 * - s1 = slope of Gibbs function at al = 0, which is the previous 885 * solution = d(Gibbs)/d(al). 886 * - s2 = slope of Gibbs function at al = 1, which is the current 887 * solution = d(Gibbs)/d(al). 888 * 889 * Only if there has been an inflection point (i.e., s1 < 0 and s2 > 0), 890 * does this code section kick in. It finds the point on the parabola 891 * where the slope is equal to zero. 892 */ 893 bool vcs_globStepDamp(); 894 895 //! Calculate the norm of a deltaGibbs free energy vector 896 /*! 897 * Positive DG for species which don't exist are ignored. 898 * 899 * @param dg Vector of local delta G's. 900 */ 901 double l2normdg(double dg[]) const; 902 903 void checkDelta1(double* const ds, double* const delTPhMoles, size_t kspec); 904 905 //! Estimate equilibrium compositions 906 /*! 907 * Algorithm covered in a section of Smith and Missen's Book. 908 * 909 * Linear programming module is based on using dbolm. 910 * 911 * @param aw aw[i[ Mole fraction work space (ne in length) 912 * @param sa sa[j] = Gram-Schmidt orthog work space (ne in length) 913 * @param sm sm[i+j*ne] = QR matrix work space (ne*ne in length) 914 * @param ss ss[j] = Gram-Schmidt orthog work space (ne in length) 915 * @param test This is a small negative number. 916 */ 917 void vcs_inest(double* const aw, double* const sa, double* const sm, 918 double* const ss, double test); 919 920 //! Calculate the status of single species phases. 921 void vcs_SSPhase(); 922 923 //! Delete memory that isn't just resizable STL containers 924 /*! 925 * This gets called by the destructor or by InitSizes(). 926 */ 927 void vcs_delete_memory(); 928 929 //! Initialize the internal counters 930 /*! 931 * Initialize the internal counters containing the subroutine call 932 * values and times spent in the subroutines. 933 * 934 * ifunc = 0 Initialize only those counters appropriate for the top of 935 * vcs_solve_TP(). 936 * = 1 Initialize all counters. 937 */ 938 void vcs_counters_init(int ifunc); 939 940 //! Create a report on the plog file containing timing and its information 941 /*! 942 * @param timing_print_lvl If 0, just report the iteration count. If larger 943 * than zero, report the timing information 944 */ 945 void vcs_TCounters_report(int timing_print_lvl = 1); 946 947 void vcs_setFlagsVolPhases(const bool upToDate, const int stateCalc); 948 949 void vcs_setFlagsVolPhase(const size_t iph, const bool upToDate, const int stateCalc); 950 951 //! Update all underlying vcs_VolPhase objects 952 /*! 953 * Update the mole numbers and the phase voltages of all phases in the vcs 954 * problem 955 * 956 * @param stateCalc Location of the update (either VCS_STATECALC_NEW or 957 * VCS_STATECALC_OLD). 958 */ 959 void vcs_updateMolNumVolPhases(const int stateCalc); 960 961 // Helper functions used internally by vcs_solve_TP 962 int solve_tp_component_calc(bool& allMinorZeroedSpecies); 963 void solve_tp_inner(size_t& iti, size_t& it1, bool& uptodate_minors, 964 bool& allMinorZeroedSpecies, int& forceComponentCalc, 965 int& stage, bool printDetails, char* ANOTE); 966 void solve_tp_equilib_check(bool& allMinorZeroedSpecies, bool& uptodate_minors, 967 bool& giveUpOnElemAbund, int& solveFail, 968 size_t& iti, size_t& it1, int maxit, 969 int& stage, bool& lec); 970 void solve_tp_elem_abund_check(size_t& iti, int& stage, bool& lec, 971 bool& giveUpOnElemAbund, 972 int& finalElemAbundAttempts, 973 int& rangeErrorFound); 974 975 // data used by vcs_solve_TP and it's helper functions 976 vector_fp m_sm; 977 vector_fp m_ss; 978 vector_fp m_sa; 979 vector_fp m_aw; 980 vector_fp m_wx; 981 982 public: 983 //! Print level for print routines 984 int m_printLvl; 985 986 MultiPhase* m_mix; 987 988 //! Print out the problem specification in all generality as it currently 989 //! exists in the VCS_SOLVE object 990 /*! 991 * @param print_lvl Parameter lvl for printing 992 * * 0 - no printing 993 * * 1 - all printing 994 */ 995 void prob_report(int print_lvl); 996 997 //! Add elements to the local element list 998 /*! 999 * This routine sorts through the elements defined in the vcs_VolPhase 1000 * object. It then adds the new elements to the VCS_SOLVE object, and creates 1001 * a global map, which is stored in the vcs_VolPhase object. Id and matching 1002 * of elements is done strictly via the element name, with case not 1003 * mattering. 1004 * 1005 * The routine also fills in the position of the element in the vcs_VolPhase 1006 * object's ElGlobalIndex field. 1007 * 1008 * @param volPhase Object containing the phase to be added. The elements in 1009 * this phase are parsed for addition to the global element list 1010 */ 1011 void addPhaseElements(vcs_VolPhase* volPhase); 1012 1013 //! This routines adds entries for the formula matrix for one species 1014 /*! 1015 * This routines adds entries for the formula matrix for this object for one 1016 * species 1017 * 1018 * This object also fills in the index filed, IndSpecies, within the 1019 * volPhase object. 1020 * 1021 * @param volPhase object containing the species 1022 * @param k Species number within the volPhase k 1023 * @param kT global Species number within this object 1024 * 1025 */ 1026 size_t addOnePhaseSpecies(vcs_VolPhase* volPhase, size_t k, size_t kT); 1027 1028 //! This routine resizes the number of elements in the VCS_SOLVE object by 1029 //! adding a new element to the end of the element list 1030 /*! 1031 * The element name is added. Formula vector entries ang element abundances 1032 * for the new element are set to zero. 1033 * 1034 * @param elNameNew New name of the element 1035 * @param elType Type of the element 1036 * @param elactive boolean indicating whether the element is active 1037 * @returns the index number of the new element 1038 */ 1039 size_t addElement(const char* elNameNew, int elType, int elactive); 1040 1041 void reportCSV(const std::string& reportFile); 1042 1043 //! Total number of species in the problems 1044 size_t m_nsp; 1045 1046 //! Number of element constraints in the problem 1047 size_t m_nelem; 1048 1049 //! Number of components calculated for the problem 1050 size_t m_numComponents; 1051 1052 //! Total number of non-component species in the problem 1053 size_t m_numRxnTot; 1054 1055 //! Current number of species in the problems. Species can be deleted if 1056 //! they aren't stable under the current conditions 1057 size_t m_numSpeciesRdc; 1058 1059 //! Current number of non-component species in the problem. Species can be 1060 //! deleted if they aren't stable under the current conditions 1061 size_t m_numRxnRdc; 1062 1063 //! Number of active species which are currently either treated as 1064 //! minor species 1065 size_t m_numRxnMinorZeroed; 1066 1067 //! Number of Phases in the problem 1068 size_t m_numPhases; 1069 1070 //! Formula matrix for the problem 1071 /*! 1072 * FormulaMatrix(kspec,j) = Number of elements, j, in the kspec species 1073 * 1074 * Both element and species indices are swapped. 1075 */ 1076 Array2D m_formulaMatrix; 1077 1078 //! Stoichiometric coefficient matrix for the reaction mechanism expressed 1079 //! in Reduced Canonical Form. 1080 /*! 1081 * This is the stoichiometric coefficient matrix for the reaction which 1082 * forms species kspec from the component species. A stoichiometric 1083 * coefficient of one is assumed for the species kspec in this mechanism. 1084 * 1085 * NOTE: kspec = irxn + m_numComponents 1086 * 1087 * `m_stoichCoeffRxnMatrix(j,irxn)` : j refers to the component number, and 1088 * irxn refers to the irxn_th non-component species. The stoichiometric 1089 * coefficients multiplied by the Formula coefficients of the component 1090 * species add up to the negative value of the number of elements in the 1091 * species kspec. 1092 * 1093 * size = nelements0 x nspecies0 1094 */ 1095 Array2D m_stoichCoeffRxnMatrix; 1096 1097 //! Absolute size of the stoichiometric coefficients 1098 /*! 1099 * scSize[irxn] = abs(Size) of the stoichiometric coefficients. These are 1100 * used to determine whether a given species should be 1101 * handled by the alt_min treatment or should be handled as a 1102 * major species. 1103 */ 1104 vector_fp m_scSize; 1105 1106 //! total size of the species 1107 /*! 1108 * This is used as a multiplier to the mole number in figuring out which 1109 * species should be components. 1110 */ 1111 vector_fp m_spSize; 1112 1113 //! Standard state chemical potentials for species K at the current 1114 //! temperature and pressure. 1115 /*! 1116 * The first NC entries are for components. The following NR entries are 1117 * for the current non-component species in the mechanism. 1118 */ 1119 vector_fp m_SSfeSpecies; 1120 1121 //! Free energy vector from the start of the current iteration 1122 /*! 1123 * The free energies are saved at the start of the current iteration. 1124 * Length = number of species 1125 */ 1126 vector_fp m_feSpecies_old; 1127 1128 //! Dimensionless new free energy for all the species in the mechanism 1129 //! at the new tentative T, P, and mole numbers. 1130 /*! 1131 * The first NC entries are for components. The following 1132 * NR entries are for the current non-component species in the mechanism. 1133 * Length = number of species 1134 */ 1135 vector_fp m_feSpecies_new; 1136 1137 //! Setting for whether to do an initial estimate 1138 /*! 1139 * Initial estimate: 1140 * 0 Do not estimate the solution at all. Use the 1141 * supplied mole numbers as is. 1142 * 1 Only do an estimate if the element abundances aren't satisfied. 1143 * -1 Force an estimate of the soln. Throw out the input mole numbers. 1144 */ 1145 int m_doEstimateEquil; 1146 1147 //! Total moles of the species 1148 /*! 1149 * Total number of moles of the kth species. 1150 * Length = Total number of species = m 1151 */ 1152 vector_fp m_molNumSpecies_old; 1153 1154 //! Specifies the species unknown type 1155 /*! 1156 * There are two types. One is the straightforward species, with the mole 1157 * number w[k], as the unknown. The second is the an interfacial voltage 1158 * where w[k] refers to the interfacial voltage in volts. 1159 * 1160 * These species types correspond to metallic electrons corresponding to 1161 * electrodes. The voltage and other interfacial conditions sets up an 1162 * interfacial current, which is set to zero in this initial treatment. 1163 * Later we may have non-zero interfacial currents. 1164 */ 1165 vector_int m_speciesUnknownType; 1166 1167 //! Change in the number of moles of phase, iphase, due to the 1168 //! noncomponent formation reaction, irxn, for species, k: 1169 /*! 1170 * m_deltaMolNumPhase(iphase,irxn) = k = nc + irxn 1171 */ 1172 Array2D m_deltaMolNumPhase; 1173 1174 //! This is 1 if the phase, iphase, participates in the formation reaction 1175 //! irxn, and zero otherwise. PhaseParticipation(iphase,irxn) 1176 Array2D m_phaseParticipation; 1177 1178 //! electric potential of the iph phase 1179 vector_fp m_phasePhi; 1180 1181 //! Tentative value of the mole number vector. It's also used to store the 1182 //! mole fraction vector. 1183 vector_fp m_molNumSpecies_new; 1184 1185 //! Delta G(irxn) for the noncomponent species in the mechanism. 1186 /*! 1187 * Computed by the subroutine deltaG. m_deltaGRxn is the free energy change 1188 * for the reaction which forms species K from the component species. This 1189 * vector has length equal to the number of noncomponent species in the 1190 * mechanism. It starts with the first current noncomponent species in the 1191 * mechanism. 1192 */ 1193 vector_fp m_deltaGRxn_new; 1194 1195 //! Last deltag[irxn] from the previous step 1196 vector_fp m_deltaGRxn_old; 1197 1198 //! Last deltag[irxn] from the previous step with additions for 1199 //! possible births of zeroed phases. 1200 vector_fp m_deltaGRxn_Deficient; 1201 1202 //! Temporary vector of Rxn DeltaG's 1203 /*! 1204 * This is used from time to time, for printing purposes 1205 */ 1206 vector_fp m_deltaGRxn_tmp; 1207 1208 //! Reaction Adjustments for each species during the current step 1209 /*! 1210 * delta Moles for each species during the current step. 1211 * Length = number of species 1212 */ 1213 vector_fp m_deltaMolNumSpecies; 1214 1215 //! Element abundances vector 1216 /*! 1217 * Vector of moles of each element actually in the solution vector. Except 1218 * for certain parts of the algorithm, this is a constant. Note other 1219 * constraint conditions are added to this vector. This is input from the 1220 * input file and is considered a constant from thereon. units = kmoles 1221 */ 1222 vector_fp m_elemAbundances; 1223 1224 //! Element abundances vector Goals 1225 /*! 1226 * Vector of moles of each element that are the goals of the simulation. 1227 * This is a constant in the problem. Note other constraint conditions are 1228 * added to this vector. This is input from the input file and is considered 1229 * a constant from thereon. units = kmoles 1230 */ 1231 vector_fp m_elemAbundancesGoal; 1232 1233 //! Total number of kmoles in all phases 1234 /*! 1235 * This number includes the inerts. 1236 * -> Don't use this except for scaling purposes 1237 */ 1238 double m_totalMolNum; 1239 1240 //! Total kmols of species in each phase 1241 /*! 1242 * This contains the total number of moles of species in each phase 1243 * 1244 * Length = number of phases 1245 */ 1246 vector_fp m_tPhaseMoles_old; 1247 1248 //! total kmols of species in each phase in the tentative soln vector 1249 /*! 1250 * This contains the total number of moles of species in each phase 1251 * in the tentative solution vector 1252 * 1253 * Length = number of phases 1254 */ 1255 vector_fp m_tPhaseMoles_new; 1256 1257 //! Temporary vector of length NPhase 1258 mutable vector_fp m_TmpPhase; 1259 1260 //! Temporary vector of length NPhase 1261 mutable vector_fp m_TmpPhase2; 1262 1263 //! Change in the total moles in each phase 1264 /*! 1265 * Length number of phases. 1266 */ 1267 vector_fp m_deltaPhaseMoles; 1268 1269 //! Temperature (Kelvin) 1270 double m_temperature; 1271 1272 //! Pressure 1273 double m_pressurePA; 1274 1275 //! Total kmoles of inert to add to each phase 1276 /*! 1277 * TPhInertMoles[iph] = Total kmoles of inert to add to each phase 1278 * length = number of phases 1279 */ 1280 vector_fp TPhInertMoles; 1281 1282 //! Tolerance requirement for major species 1283 double m_tolmaj; 1284 1285 //! Tolerance requirements for minor species 1286 double m_tolmin; 1287 1288 //! Below this, major species aren't refined any more 1289 double m_tolmaj2; 1290 1291 //! Below this, minor species aren't refined any more 1292 double m_tolmin2; 1293 1294 //! Index vector that keeps track of the species vector rearrangement 1295 /*! 1296 * At the end of each run, the species vector and associated data gets put 1297 * back in the original order. 1298 * 1299 * Example 1300 * 1301 * k = m_speciesMapIndex[kspec] 1302 * 1303 * kspec = current order in the vcs_solve object 1304 * k = original order in the MultiPhase object 1305 */ 1306 std::vector<size_t> m_speciesMapIndex; 1307 1308 //! Index that keeps track of the index of the species within the local 1309 //! phase 1310 /*! 1311 * This returns the local index of the species within the phase. Its 1312 * argument is the global species index within the VCS problem. 1313 * 1314 * k = m_speciesLocalPhaseIndex[kspec] 1315 * 1316 * k varies between 0 and the nSpecies in the phase 1317 * 1318 * Length = number of species 1319 */ 1320 std::vector<size_t> m_speciesLocalPhaseIndex; 1321 1322 //! Index vector that keeps track of the rearrangement of the elements 1323 /*! 1324 * At the end of each run, the element vector and associated data gets 1325 * put back in the original order. 1326 * 1327 * Example 1328 * 1329 * e = m_elementMapIndex[eNum] 1330 * eNum = current order in the vcs_solve object 1331 * e = original order in the MultiPhase object 1332 */ 1333 std::vector<size_t> m_elementMapIndex; 1334 1335 //! Mapping between the species index for noncomponent species and the 1336 //! full species index. 1337 /*! 1338 * ir[irxn] = Mapping between the reaction index for noncomponent 1339 * formation reaction of a species and the full species 1340 * index. 1341 * 1342 * Initially set to a value of K = NC + I This vector has length equal to 1343 * number of noncomponent species in the mechanism. It starts with the 1344 * first current noncomponent species in the mechanism. kspec = ir[irxn] 1345 */ 1346 std::vector<size_t> m_indexRxnToSpecies; 1347 1348 //! Major -Minor status vector for the species in the problem 1349 /*! 1350 * The index for this is species. The reaction that this is referring to 1351 * is `kspec = irxn + m_numComponents`. For possible values and their 1352 * meanings, see vcs_evaluate_speciesType(). 1353 */ 1354 vector_int m_speciesStatus; 1355 1356 //! Mapping from the species number to the phase number 1357 std::vector<size_t> m_phaseID; 1358 1359 //! Boolean indicating whether a species belongs to a single-species phase 1360 // vector<bool> can't be used here because it doesn't work with std::swap 1361 std::vector<char> m_SSPhase; 1362 1363 //! Species string name for the kth species 1364 std::vector<std::string> m_speciesName; 1365 1366 //! Vector of strings containing the element names 1367 /*! 1368 * ElName[j] = String containing element names 1369 */ 1370 std::vector<std::string> m_elementName; 1371 1372 //! Type of the element constraint 1373 /*! 1374 * * 0 - #VCS_ELEM_TYPE_ABSPOS Normal element that is positive or zero in 1375 * all species. 1376 * * 1 - #VCS_ELEM_TYPE_ELECTRONCHARGE element dof that corresponds to 1377 * the electronic charge DOF. 1378 * * 2 - #VCS_ELEM_TYPE_CHARGENEUTRALITY element dof that corresponds to 1379 * a required charge neutrality constraint on the phase. The element 1380 * abundance is always exactly zero. 1381 * * 3 - #VCS_ELEM_TYPE_OTHERCONSTRAINT Other constraint which may mean 1382 * that a species has neg 0 or pos value of that constraint (other than 1383 * charge) 1384 */ 1385 vector_int m_elType; 1386 1387 //! Specifies whether an element constraint is active 1388 /*! 1389 * The default is true. Length = nelements 1390 */ 1391 vector_int m_elementActive; 1392 1393 //! Array of Phase Structures. Length = number of phases. 1394 std::vector<std::unique_ptr<vcs_VolPhase>> m_VolPhaseList; 1395 1396 //! specifies the activity convention of the phase containing the species 1397 /*! 1398 * * 0 = molar based 1399 * * 1 = molality based 1400 * 1401 * length = number of species 1402 */ 1403 vector_int m_actConventionSpecies; 1404 1405 //! specifies the activity convention of the phase. 1406 /*! 1407 * * 0 = molar based 1408 * * 1 = molality based 1409 * 1410 * length = number of phases 1411 */ 1412 vector_int m_phaseActConvention; 1413 1414 //! specifies the ln(Mnaught) used to calculate the chemical potentials 1415 /*! 1416 * For molar based activity conventions this will be equal to 0.0. 1417 * length = number of species. 1418 */ 1419 vector_fp m_lnMnaughtSpecies; 1420 1421 //! Molar-based Activity Coefficients for Species. 1422 //! Length = number of species 1423 vector_fp m_actCoeffSpecies_new; 1424 1425 //! Molar-based Activity Coefficients for Species based on old mole numbers 1426 /*! 1427 * These activity coefficients are based on the m_molNumSpecies_old 1428 * values Molar based activity coefficients. Length = number of species 1429 */ 1430 vector_fp m_actCoeffSpecies_old; 1431 1432 //! Change in the log of the activity coefficient with respect to the mole 1433 //! number multiplied by the phase mole number 1434 /*! 1435 * size = nspecies x nspecies 1436 * 1437 * This is a temporary array that gets regenerated every time it's needed. 1438 * It is not swapped wrt species. 1439 */ 1440 Array2D m_np_dLnActCoeffdMolNum; 1441 1442 //! Molecular weight of each species 1443 /*! 1444 * units = kg/kmol. length = number of species. 1445 * 1446 * note: this is a candidate for removal. I don't think we use it. 1447 */ 1448 vector_fp m_wtSpecies; 1449 1450 //! Charge of each species. Length = number of species. 1451 vector_fp m_chargeSpecies; 1452 1453 std::vector<std::vector<size_t> > phasePopProblemLists_; 1454 1455 //! Vector of pointers to thermo structures which identify the model 1456 //! and parameters for evaluating the thermodynamic functions for that 1457 //! particular species. 1458 /*! 1459 * SpeciesThermo[k] pointer to the thermo information for the kth species 1460 */ 1461 std::vector<std::unique_ptr<VCS_SPECIES_THERMO>> m_speciesThermoList; 1462 1463 //! Choice of Hessians 1464 /*! 1465 * If this is true, then we will use a better approximation to the Hessian 1466 * based on Jacobian of the ln(ActCoeff) with respect to mole numbers 1467 */ 1468 int m_useActCoeffJac; 1469 1470 //! Total volume of all phases. Units are m^3 1471 double m_totalVol; 1472 1473 //! Partial molar volumes of the species 1474 /*! 1475 * units = mks (m^3/kmol) 1476 * Length = number of species 1477 */ 1478 vector_fp m_PMVolumeSpecies; 1479 1480 //! dimensionless value of Faraday's constant, F / RT (1/volt) 1481 double m_Faraday_dim; 1482 1483 //! Timing and iteration counters for the vcs object 1484 VCS_COUNTERS* m_VCount; 1485 1486 //! Debug printing lvl 1487 /*! 1488 * Levels correspond to the following guidelines 1489 * * 0 No printing at all 1490 * * 1 Serious warnings or fatal errors get one line 1491 * * 2 one line per each successful vcs package call 1492 * * 3 one line per every successful solve_TP calculation 1493 * * 4 one line for every successful operation -> solve_TP gets a summary report 1494 * * 5 each iteration in solve_TP gets a report with one line per species 1495 * * 6 Each decision in solve_TP gets a line per species in addition to 4 1496 * * 10 Additionally Hessian matrix is printed out 1497 */ 1498 int m_debug_print_lvl; 1499 1500 //! printing level of timing information 1501 /*! 1502 * * 1 allowing printing of timing 1503 * * 0 do not allow printing of timing -> everything is printed as a NA. 1504 */ 1505 int m_timing_print_lvl; 1506 1507 //! Disable printing of timing information. Used to generate consistent 1508 //! output for tests. 1509 static void disableTiming(); 1510 1511 friend class vcs_phaseStabilitySolve; 1512 }; 1513 1514 } 1515 #endif 1516