1 /** 2 * @file MultiPhase.h 3 * Headers for the \link Cantera::MultiPhase MultiPhase\endlink 4 * object that is used to set up multiphase equilibrium problems (see \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 CT_MULTIPHASE_H 11 #define CT_MULTIPHASE_H 12 13 #include "cantera/numerics/DenseMatrix.h" 14 15 namespace Cantera 16 { 17 18 class ThermoPhase; 19 20 //! @defgroup equilfunctions 21 22 //! A class for multiphase mixtures. The mixture can contain any 23 //! number of phases of any type. 24 /*! 25 * This object is the basic tool used by Cantera for use in Multiphase 26 * equilibrium calculations. 27 * 28 * It is a container for a set of phases. Each phase has a given number of 29 * kmoles. Therefore, MultiPhase may be considered an "extrinsic" 30 * thermodynamic object, in contrast to the ThermoPhase object, which is an 31 * "intrinsic" thermodynamic object. 32 * 33 * MultiPhase may be considered to be "upstream" of the ThermoPhase objects in 34 * the sense that setting a property within MultiPhase, such as temperature, 35 * pressure, or species mole number, affects the underlying ThermoPhase 36 * object, but not the other way around. 37 * 38 * All phases have the same temperature and pressure, and a specified number 39 * of moles for each phase. The phases do not need to have the same elements. 40 * For example, a mixture might consist of a gaseous phase with elements (H, 41 * C, O, N), a solid carbon phase containing only element C, etc. A master 42 * element set will be constructed for the mixture that is the intersection of 43 * the elements of each phase. 44 * 45 * Below, reference is made to global species and global elements. These refer 46 * to the collective species and elements encompassing all of the phases 47 * tracked by the object. 48 * 49 * The global element list kept by this object is an intersection of the 50 * element lists of all the phases that comprise the MultiPhase. 51 * 52 * The global species list kept by this object is a concatenated list of all 53 * of the species in all the phases that comprise the MultiPhase. The ordering 54 * of species is contiguous with respect to the phase id. 55 * 56 * @ingroup equilfunctions 57 */ 58 class MultiPhase 59 { 60 public: 61 //! Constructor. 62 /*! 63 * The constructor takes no arguments, since phases are added using 64 * method addPhase(). 65 */ 66 MultiPhase(); 67 68 //! Destructor. Does nothing. Class MultiPhase does not take "ownership" 69 //! (i.e. responsibility for destroying) the phase objects. ~MultiPhase()70 virtual ~MultiPhase() {} 71 72 //! Add a vector of phases to the mixture 73 /*! 74 * See the single addPhases command. This just does a bunch of phases 75 * at one time 76 * @param phases Vector of pointers to phases 77 * @param phaseMoles Vector of mole numbers in each phase (kmol) 78 */ 79 void addPhases(std::vector<ThermoPhase*>& phases, const vector_fp& phaseMoles); 80 81 //! Add all phases present in 'mix' to this mixture. 82 /*! 83 * @param mix Add all of the phases in another MultiPhase 84 * object to the current object. 85 */ 86 void addPhases(MultiPhase& mix); 87 88 //! Add a phase to the mixture. 89 /*! 90 * This function must be called before the init() function is called, 91 * which serves to freeze the MultiPhase. 92 * 93 * @param p pointer to the phase object 94 * @param moles total number of moles of all species in this phase 95 */ 96 void addPhase(ThermoPhase* p, doublereal moles); 97 98 //! Number of elements. nElements()99 size_t nElements() const { 100 return m_nel; 101 } 102 103 //! Check that the specified element index is in range. 104 //! Throws an exception if m is greater than nElements()-1 105 void checkElementIndex(size_t m) const; 106 107 //! Check that an array size is at least nElements(). 108 //! Throws an exception if mm is less than nElements(). Used before calls 109 //! which take an array pointer. 110 void checkElementArraySize(size_t mm) const; 111 112 //! Returns the name of the global element *m*. 113 /*! 114 * @param m index of the global element 115 */ 116 std::string elementName(size_t m) const; 117 118 //! Returns the index of the element with name \a name. 119 /*! 120 * @param name String name of the global element 121 */ 122 size_t elementIndex(const std::string& name) const; 123 124 //! Number of species, summed over all phases. nSpecies()125 size_t nSpecies() const { 126 return m_nsp; 127 } 128 129 //! Check that the specified species index is in range. 130 //! Throws an exception if k is greater than nSpecies()-1 131 void checkSpeciesIndex(size_t k) const; 132 133 //! Check that an array size is at least nSpecies(). 134 //! Throws an exception if kk is less than nSpecies(). Used before calls 135 //! which take an array pointer. 136 void checkSpeciesArraySize(size_t kk) const; 137 138 //! Name of species with global index \a kGlob 139 /*! 140 * @param kGlob global species index 141 */ 142 std::string speciesName(const size_t kGlob) const; 143 144 //! Returns the Number of atoms of global element \a mGlob in 145 //! global species \a kGlob. 146 /*! 147 * @param kGlob global species index 148 * @param mGlob global element index 149 * @returns the number of atoms. 150 */ 151 doublereal nAtoms(const size_t kGlob, const size_t mGlob) const; 152 153 /// Returns the global Species mole fractions. 154 /*! 155 * Write the array of species mole 156 * fractions into array \c x. The mole fractions are 157 * normalized to sum to one in each phase. 158 * 159 * @param x vector of mole fractions. Length = number of global species. 160 */ 161 void getMoleFractions(doublereal* const x) const; 162 163 //! Process phases and build atomic composition array. 164 /*! 165 * This method must be called after all phases are added, before doing 166 * anything else with the mixture. After init() has been called, no more 167 * phases may be added. 168 */ 169 void init(); 170 171 //! Returns the name of the n'th phase 172 /*! 173 * @param iph phase Index 174 */ 175 std::string phaseName(const size_t iph) const; 176 177 //! Returns the index, given the phase name 178 /*! 179 * @param pName Name of the phase 180 * @returns the index. A value of -1 means the phase isn't in the object. 181 */ 182 int phaseIndex(const std::string& pName) const; 183 184 //! Return the number of moles in phase n. 185 /*! 186 * @param n Index of the phase. 187 */ 188 doublereal phaseMoles(const size_t n) const; 189 190 //! Set the number of moles of phase with index n. 191 /*! 192 * @param n Index of the phase 193 * @param moles Number of moles in the phase (kmol) 194 */ 195 void setPhaseMoles(const size_t n, const doublereal moles); 196 197 /// Return a reference to phase n. 198 /*! 199 * The state of phase n is also updated to match the state stored locally 200 * in the mixture object. 201 * 202 * @param n Phase Index 203 * @return Reference to the ThermoPhase object for the phase 204 */ 205 ThermoPhase& phase(size_t n); 206 207 //! Check that the specified phase index is in range 208 //! Throws an exception if m is greater than nPhases() 209 void checkPhaseIndex(size_t m) const; 210 211 //! Check that an array size is at least nPhases() 212 //! Throws an exception if mm is less than nPhases(). Used before calls 213 //! which take an array pointer. 214 void checkPhaseArraySize(size_t mm) const; 215 216 //! Returns the moles of global species \c k. units = kmol 217 /*! 218 * @param kGlob Global species index k 219 */ 220 doublereal speciesMoles(size_t kGlob) const; 221 222 //! Return the global index of the species belonging to phase number \c p 223 //! with local index \c k within the phase. 224 /*! 225 * @param k local index of the species within the phase 226 * @param p index of the phase 227 */ speciesIndex(size_t k,size_t p)228 size_t speciesIndex(size_t k, size_t p) const { 229 return m_spstart[p] + k; 230 } 231 232 //! Return the global index of the species belonging to phase name \c phaseName 233 //! with species name \c speciesName 234 /*! 235 * @param speciesName Species Name 236 * @param phaseName Phase Name 237 * 238 * @returns the global index 239 * 240 * If the species or phase name is not recognized, this routine throws a 241 * CanteraError. 242 */ 243 size_t speciesIndex(const std::string& speciesName, const std::string& phaseName); 244 245 /// Minimum temperature for which all solution phases have valid thermo 246 /// data. Stoichiometric phases are not considered, since they may have 247 /// thermo data only valid for conditions for which they are stable. minTemp()248 doublereal minTemp() const { 249 return m_Tmin; 250 } 251 252 /// Maximum temperature for which all solution phases have valid thermo 253 /// data. Stoichiometric phases are not considered, since they may have 254 /// thermo data only valid for conditions for which they are stable. maxTemp()255 doublereal maxTemp() const { 256 return m_Tmax; 257 } 258 259 //! Total charge summed over all phases (Coulombs). 260 doublereal charge() const; 261 262 /// Charge (Coulombs) of phase with index \a p. 263 /*! 264 * The net charge is computed as \f[ Q_p = N_p \sum_k F z_k X_k \f] 265 * where the sum runs only over species in phase \a p. 266 * @param p index of the phase for which the charge is desired. 267 */ 268 doublereal phaseCharge(size_t p) const; 269 270 //! Total moles of global element \a m, summed over all phases. 271 /*! 272 * @param m Index of the global element 273 */ 274 doublereal elementMoles(size_t m) const; 275 276 //! Returns a vector of Chemical potentials. 277 /*! 278 * Write into array \a mu the chemical potentials of all species 279 * [J/kmol]. The chemical potentials are related to the activities by 280 * 281 * \f$ 282 * \mu_k = \mu_k^0(T, P) + RT \ln a_k. 283 * \f$. 284 * 285 * @param mu Chemical potential vector. Length = num global species. Units 286 * = J/kmol. 287 */ 288 void getChemPotentials(doublereal* mu) const; 289 290 /// Returns a vector of Valid chemical potentials. 291 /*! 292 * Write into array \a mu the chemical potentials of all species with 293 * thermo data valid for the current temperature [J/kmol]. For other 294 * species, set the chemical potential to the value \a not_mu. If \a 295 * standard is set to true, then the values returned are standard chemical 296 * potentials. 297 * 298 * This method is designed for use in computing chemical equilibrium by 299 * Gibbs minimization. For solution phases (more than one species), this 300 * does the same thing as getChemPotentials. But for stoichiometric 301 * phases, this writes into array \a mu the user-specified value \a not_mu 302 * instead of the chemical potential if the temperature is outside the 303 * range for which the thermo data for the one species in the phase are 304 * valid. The need for this arises since many condensed phases have thermo 305 * data fit only for the temperature range for which they are stable. For 306 * example, in the NASA database, the fits for H2O(s) are only done up to 307 * 0 C, the fits for H2O(L) are only done from 0 C to 100 C, etc. Using 308 * the polynomial fits outside the range for which the fits were done can 309 * result in spurious chemical potentials, and can lead to condensed 310 * phases appearing when in fact they should be absent. 311 * 312 * By setting \a not_mu to a large positive value, it is possible to force 313 * routines which seek to minimize the Gibbs free energy of the mixture to 314 * zero out any phases outside the temperature range for which their 315 * thermo data are valid. 316 * 317 * @param not_mu Value of the chemical potential to set species in phases, 318 * for which the thermo data is not valid 319 * @param mu Vector of chemical potentials. length = Global species, 320 * units = J kmol-1 321 * @param standard If this method is called with \a standard set to true, 322 * then the composition-independent standard chemical 323 * potentials are returned instead of the composition- 324 * dependent chemical potentials. 325 */ 326 void getValidChemPotentials(doublereal not_mu, doublereal* mu, 327 bool standard = false) const; 328 329 //! Temperature [K]. temperature()330 doublereal temperature() const { 331 return m_temp; 332 } 333 334 //! Equilibrate a MultiPhase object 335 /*! 336 * Set this mixture to chemical equilibrium by calling one of Cantera's 337 * equilibrium solvers. The XY parameter indicates what two thermodynamic 338 * quantities are to be held constant during the equilibration process. 339 * 340 * @param XY String representation of what two properties are being 341 * held constant 342 * @param solver Name of the solver to be used to equilibrate the phase. 343 * If solver = 'vcs', the vcs_MultiPhaseEquil solver will be used. If 344 * solver = 'gibbs', the MultiPhaseEquil solver will be used. If solver 345 * = 'auto', the 'vcs' solver will be tried first, followed by the 346 * 'gibbs' solver if the first one fails. 347 * @param rtol Relative tolerance 348 * @param max_steps Maximum number of steps to take to find the solution 349 * @param max_iter The maximum number of outer temperature or pressure 350 * iterations to take when T and/or P is not held fixed. 351 * @param estimate_equil integer indicating whether the solver should 352 * estimate its own initial condition. If 0, the initial mole fraction 353 * vector in the ThermoPhase object is used as the initial condition. 354 * If 1, the initial mole fraction vector is used if the element 355 * abundances are satisfied. If -1, the initial mole fraction vector is 356 * thrown out, and an estimate is formulated. 357 * @param log_level loglevel Controls amount of diagnostic output. 358 * log_level=0 suppresses diagnostics, and increasingly-verbose 359 * messages are written as loglevel increases. 360 * 361 * @ingroup equilfunctions 362 */ 363 void equilibrate(const std::string& XY, const std::string& solver="auto", 364 double rtol=1e-9, int max_steps=50000, int max_iter=100, 365 int estimate_equil=0, int log_level=0); 366 367 /// Set the temperature [K]. 368 /*! 369 * @param T value of the temperature (Kelvin) 370 */ 371 void setTemperature(const doublereal T); 372 373 //! Set the state of the underlying ThermoPhase objects in one call 374 /*! 375 * @param T Temperature of the system (kelvin) 376 * @param Pres pressure of the system (pascal) 377 */ 378 void setState_TP(const doublereal T, const doublereal Pres); 379 380 //! Set the state of the underlying ThermoPhase objects in one call 381 /*! 382 * @param T Temperature of the system (kelvin) 383 * @param Pres pressure of the system (pascal) 384 * @param Moles Vector of mole numbers of all the species in all the phases 385 * (kmol) 386 */ 387 void setState_TPMoles(const doublereal T, const doublereal Pres, const doublereal* Moles); 388 389 /// Pressure [Pa]. pressure()390 doublereal pressure() const { 391 return m_press; 392 } 393 394 /// The total mixture volume [m^3]. 395 /*! 396 * Returns the cumulative sum of the volumes of all the phases in the 397 * mixture. 398 */ 399 doublereal volume() const; 400 401 //! Set the pressure [Pa]. 402 /*! 403 * @param P Set the pressure in the MultiPhase object (Pa) 404 */ setPressure(doublereal P)405 void setPressure(doublereal P) { 406 m_press = P; 407 updatePhases(); 408 } 409 410 //! The enthalpy of the mixture [J]. 411 doublereal enthalpy() const; 412 413 //! The internal energy of the mixture [J]. 414 doublereal IntEnergy() const; 415 416 //! The entropy of the mixture [J/K]. 417 doublereal entropy() const; 418 419 //! The Gibbs function of the mixture [J]. 420 doublereal gibbs() const; 421 422 //! Heat capacity at constant pressure [J/K]. Note that this does not 423 //! account for changes in composition of the mixture with temperature. 424 doublereal cp() const; 425 426 //! Number of phases. nPhases()427 size_t nPhases() const { 428 return m_phase.size(); 429 } 430 431 //! Return true is species \a kGlob is a species in a multicomponent 432 //! solution phase. 433 /*! 434 * @param kGlob index of the global species 435 */ 436 bool solutionSpecies(size_t kGlob) const; 437 438 //! Returns the phase index of the Kth "global" species 439 /*! 440 * @param kGlob Global species index. 441 * @returns the index of the owning phase. 442 */ 443 size_t speciesPhaseIndex(const size_t kGlob) const; 444 445 //! Returns the mole fraction of global species k 446 /*! 447 * @param kGlob Index of the global species. 448 */ 449 doublereal moleFraction(const size_t kGlob) const; 450 451 //! Set the Mole fractions of the nth phase 452 /*! 453 * This function sets the mole fractions of the nth phase. Note, the mole 454 * number of the phase stays constant 455 * 456 * @param n index of the phase 457 * @param x Vector of input mole fractions. 458 */ 459 void setPhaseMoleFractions(const size_t n, const doublereal* const x); 460 461 //! Set the number of moles of species in the mixture 462 /*! 463 * @param xMap CompositionMap of the species with nonzero mole numbers. 464 * Mole numbers that are less than or equal to zero will be 465 * set to zero. units = kmol. 466 */ 467 void setMolesByName(const compositionMap& xMap); 468 469 //! Set the moles via a string containing their names. 470 /*! 471 * The string x is in the form of a composition map. Species which are not 472 * listed are set to zero. 473 * 474 * @param x string x in the form of a composition map 475 * where values are the moles of the species. 476 */ 477 void setMolesByName(const std::string& x); 478 479 //! Get the mole numbers of all species in the multiphase object 480 /*! 481 * @param[out] molNum Vector of doubles of length nSpecies containing the 482 * global mole numbers (kmol). 483 */ 484 void getMoles(doublereal* molNum) const; 485 486 //! Sets all of the global species mole numbers 487 /*! 488 * The state of each phase object is also updated to have the specified 489 * composition and the mixture temperature and pressure. 490 * 491 * @param n Vector of doubles of length nSpecies containing the global 492 * mole numbers (kmol). 493 */ 494 void setMoles(const doublereal* n); 495 496 //! Adds moles of a certain species to the mixture 497 /*! 498 * @param indexS Index of the species in the MultiPhase object 499 * @param addedMoles Value of the moles that are added to the species. 500 */ 501 void addSpeciesMoles(const int indexS, const doublereal addedMoles); 502 503 //! Retrieves a vector of element abundances 504 /*! 505 * @param elemAbundances Vector of element abundances 506 * Length = number of elements in the MultiPhase object. 507 * Index is the global element index. Units is in kmol. 508 */ 509 void getElemAbundances(doublereal* elemAbundances) const; 510 511 //! Return true if the phase \a p has valid thermo data for the current 512 //! temperature. 513 /*! 514 * @param p Index of the phase. 515 */ 516 bool tempOK(size_t p) const; 517 518 // These methods are meant for internal use. 519 520 //! Update the locally-stored composition within this object to match the 521 //! current compositions of the phase objects. 522 /*! 523 * Query the underlying ThermoPhase objects for their mole fractions and 524 * fill in the mole fraction vector of this current object. Adjust element 525 * compositions within this object to match. 526 * 527 * This is an upload operation in the sense that we are taking downstream 528 * information (ThermoPhase object info) and applying it to an upstream 529 * object (MultiPhase object). 530 */ 531 void uploadMoleFractionsFromPhases(); 532 533 //! Set the states of the phase objects to the locally-stored 534 //! state within this MultiPhase object. 535 /*! 536 * This method sets each phase to the mixture temperature and pressure, 537 * and sets the phase mole fractions based on the mixture mole numbers. 538 * 539 * This is an download operation in the sense that we are taking upstream 540 * object information (MultiPhase object) and applying it to downstream 541 * objects (ThermoPhase object information) 542 * 543 * Therefore, the term, "update", is appropriate for a downstream operation. 544 */ 545 void updatePhases() const; 546 547 private: 548 //! Calculate the element abundance vector 549 void calcElemAbundances() const; 550 551 //! Set the mixture to a state of chemical equilibrium using the 552 //! MultiPhaseEquil solver. 553 /*! 554 * @param XY Integer flag specifying properties to hold fixed. 555 * @param err Error tolerance for \f$\Delta \mu/RT \f$ for all reactions. 556 * Also used as the relative error tolerance for the outer loop. 557 * @param maxsteps Maximum number of steps to take in solving the fixed TP 558 * problem. 559 * @param maxiter Maximum number of "outer" iterations for problems holding 560 * fixed something other than (T,P). 561 * @param loglevel Level of diagnostic output 562 */ 563 double equilibrate_MultiPhaseEquil(int XY, doublereal err, int maxsteps, 564 int maxiter, int loglevel); 565 566 //! Vector of the number of moles in each phase. 567 /*! 568 * Length = m_np, number of phases. 569 */ 570 vector_fp m_moles; 571 572 //! Vector of the ThermoPhase pointers. 573 std::vector<ThermoPhase*> m_phase; 574 575 //! Global Stoichiometric Coefficient array 576 /*! 577 * This is a two dimensional array m_atoms(m, k). The first index is the 578 * global element index. The second index, k, is the global species index. 579 * The value is the number of atoms of type m in species k. 580 */ 581 DenseMatrix m_atoms; 582 583 //! Locally stored vector of mole fractions of all species comprising the 584 //! MultiPhase object. 585 vector_fp m_moleFractions; 586 587 //! Mapping between the global species number and the phase ID 588 /*! 589 * m_spphase[kGlobal] = iPhase 590 * Length = number of global species 591 */ 592 std::vector<size_t> m_spphase; 593 594 //! Vector of ints containing of first species index in the global list of 595 //! species for each phase 596 /*! 597 * kfirst = m_spstart[ip], kfirst is the index of the first species in 598 * the ip'th phase. 599 */ 600 std::vector<size_t> m_spstart; 601 602 //! String names of the global elements. This has a length equal to the 603 //! number of global elements. 604 std::vector<std::string> m_enames; 605 606 //! Atomic number of each global element. 607 vector_int m_atomicNumber; 608 609 //! Vector of species names in the problem. Vector is over all species 610 //! defined in the object, the global species index. 611 std::vector<std::string> m_snames; 612 613 //! Returns the global element index, given the element string name 614 /*! 615 * -> used in the construction. However, wonder if it needs to be global. 616 */ 617 std::map<std::string, size_t> m_enamemap; 618 619 //! Current value of the temperature (kelvin) 620 doublereal m_temp; 621 622 //! Current value of the pressure (Pa) 623 doublereal m_press; 624 625 //! Number of distinct elements in all of the phases 626 size_t m_nel; 627 628 //! Number of distinct species in all of the phases 629 size_t m_nsp; 630 631 //! True if the init() routine has been called, and the MultiPhase frozen 632 bool m_init; 633 634 //! Global ID of the element corresponding to the electronic charge. If 635 //! there is none, then this is equal to -1 636 size_t m_eloc; 637 638 //! Vector of bools indicating whether temperatures are ok for phases. 639 /*! 640 * If the current temperature is outside the range of valid temperatures 641 * for the phase thermodynamics, the phase flag is set to false. 642 */ 643 mutable std::vector<bool> m_temp_OK; 644 645 //! Minimum temperature for which thermo parameterizations are valid. 646 //! Stoichiometric phases are ignored in this determination. units Kelvin 647 doublereal m_Tmin; 648 649 //! Minimum temperature for which thermo parameterizations are valid. 650 //! Stoichiometric phases are ignored in this determination. units Kelvin 651 doublereal m_Tmax; 652 653 //! Vector of element abundances 654 /*! 655 * m_elemAbundances[mGlobal] = kmol of element mGlobal summed over all 656 * species in all phases. 657 */ 658 mutable vector_fp m_elemAbundances; 659 }; 660 661 //! Function to output a MultiPhase description to a stream 662 /*! 663 * Writes out a description of the contents of each phase of the 664 * MultiPhase using the report function. 665 * 666 * @param s ostream 667 * @param x Reference to a MultiPhase 668 * @returns a reference to the ostream 669 */ 670 std::ostream& operator<<(std::ostream& s, MultiPhase& x); 671 672 //! Choose the optimum basis of species for the equilibrium calculations. 673 /*! 674 * This is done by choosing the species with the largest mole fraction not 675 * currently a linear combination of the previous components. Then, calculate 676 * the stoichiometric coefficient matrix for that basis. 677 * 678 * Calculates the identity of the component species in the mechanism. Rearranges 679 * the solution data to put the component data at the front of the species list. 680 * 681 * Then, calculates SC(J,I) the formation reactions for all noncomponent 682 * species in the mechanism. 683 * 684 * @param[in] mphase Pointer to the multiphase object. Contains the species 685 * mole fractions, which are used to pick the current optimal species 686 * component basis. 687 * @param[in] orderVectorElements Order vector for the elements. The element 688 * rows in the formula matrix are rearranged according to this vector. 689 * @param[in] orderVectorSpecies Order vector for the species. The species are 690 * rearranged according to this formula. The first nCompoments of this 691 * vector contain the calculated species components on exit. 692 * @param[in] doFormRxn If true, the routine calculates the formation 693 * reaction matrix based on the calculated component species. If 694 * false, this step is skipped. 695 * @param[out] usedZeroedSpecies = If true, then a species with a zero 696 * concentration was used as a component. The problem may be converged. 697 * @param[out] formRxnMatrix 698 * @return The number of components. 699 * 700 * @ingroup equilfunctions 701 */ 702 size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, 703 MultiPhase* mphase, std::vector<size_t>& orderVectorSpecies, 704 std::vector<size_t>& orderVectorElements, 705 vector_fp& formRxnMatrix); 706 707 //! Handles the potential rearrangement of the constraint equations 708 //! represented by the Formula Matrix. 709 /*! 710 * Rearrangement is only necessary when the number of components is less 711 * than the number of elements. For this case, some constraints can never 712 * be satisfied exactly, because the range space represented by the Formula 713 * Matrix of the components can't span the extra space. These constraints, 714 * which are out of the range space of the component Formula matrix 715 * entries, are migrated to the back of the Formula matrix. 716 * 717 * A prototypical example is an extra element column in FormulaMatrix[], which 718 * is identically zero. For example, let's say that argon is has an element 719 * column in FormulaMatrix[], but no species in the mechanism actually 720 * contains argon. Then, nc < ne. Unless the entry for desired element 721 * abundance vector for Ar is zero, then this element abundance constraint can 722 * never be satisfied. The constraint vector is not in the range space of the 723 * formula matrix. 724 * 725 * Also, without perturbation of FormulaMatrix[], BasisOptimize[] would 726 * produce a zero pivot because the matrix would be singular (unless the argon 727 * element column was already the last column of FormulaMatrix[]. 728 * 729 * This routine borrows heavily from BasisOptimize algorithm. It finds nc 730 * constraints which span the range space of the Component Formula matrix, and 731 * assigns them as the first nc components in the formula matrix. This 732 * guarantees that BasisOptimize has a nonsingular matrix to invert. 733 * 734 * @param[in] nComponents Number of components calculated previously. 735 * @param[in] elementAbundances Current value of the element abundances 736 * @param[in] mphase Input pointer to a MultiPhase object 737 * @param[in] orderVectorSpecies input vector containing the ordering of the 738 * global species in mphase. This is used to extract the component 739 * basis of the mphase object. 740 * @param[out] orderVectorElements Output vector containing the order of the 741 * elements that is necessary for calculation of the formula matrix. 742 * 743 * @ingroup equilfunctions 744 */ 745 void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances, 746 MultiPhase* mphase, 747 std::vector<size_t>& orderVectorSpecies, 748 std::vector<size_t>& orderVectorElements); 749 750 //! External int that is used to turn on debug printing for the 751 //! BasisOptimze program. 752 /*! 753 * Set this to 1 if you want debug printing from BasisOptimize. 754 */ 755 extern int BasisOptimize_print_lvl; 756 } 757 758 #endif 759