00001 /** 00002 * @file vcs_MultiPhaseEquil.h 00003 * Interface class for the vcsnonlinear solver 00004 */ 00005 00006 /* 00007 * $Author: hkmoffa $ 00008 * $Revision: 368 $ 00009 * $Date: 2010-01-03 19:46:26 -0500 (Sun, 03 Jan 2010) $ 00010 */ 00011 00012 00013 #ifndef VCS_MULTIPHASEEQUIL_H 00014 #define VCS_MULTIPHASEEQUIL_H 00015 00016 #include "ct_defs.h" 00017 #include "MultiPhase.h" 00018 00019 00020 namespace Cantera { 00021 00022 00023 //! Set a single-phase chemical solution to chemical equilibrium. 00024 /*! 00025 * The function uses the element abundance vector that is 00026 * currently consistent with the composition within the phase 00027 * itself. Two other thermodynamic quantities, determined by the 00028 * XY string, are held constant during the equilibration. 00029 * This is a convenience function that uses one or the other of 00030 * the two chemical equilibrium solvers. 00031 * 00032 * @param s The object to set to an equilibrium state 00033 * 00034 * @param XY An integer specifying the two properties to be held 00035 * constant. 00036 * 00037 * @param estimateEquil integer indicating whether the solver 00038 * should estimate its own initial condition. 00039 * If 0, the initial mole fraction vector 00040 * in the %ThermoPhase object is used as the 00041 * initial condition. 00042 * If 1, the initial mole fraction vector 00043 * is used if the element abundances are 00044 * satisfied. 00045 * if -1, the initial mole fraction vector 00046 * is thrown out, and an estimate is 00047 * formulated. 00048 * 00049 * @param printLvl Determines the amount of printing that 00050 * gets sent to stdout from the vcs package 00051 * (Note, you may have to compile with debug 00052 * flags to get some printing). 00053 * 00054 * @param solver The equilibrium solver to use. If solver = 0, 00055 * the ChemEquil solver will be used, and if 00056 * solver = 1, the vcs_MultiPhaseEquil solver will 00057 * be used (slower than ChemEquil, 00058 * but more stable). If solver < 0 (default, then 00059 * ChemEquil will be tried first, and if it fails 00060 * vcs_MultiPhaseEquil will be tried. 00061 * 00062 * @param rtol Relative tolerance of the solve. Defaults to 00063 * 1.0E-9. 00064 * 00065 * @param maxsteps The maximum number of steps to take to find 00066 * the solution. 00067 * 00068 * @param maxiter For the MultiPhaseEquil solver only, this is 00069 * the maximum number of outer temperature or 00070 * pressure iterations to take when T and/or P is 00071 * not held fixed. 00072 * 00073 * @param loglevel Controls amount of diagnostic output. loglevel 00074 * = 0 suppresses diagnostics, and increasingly-verbose 00075 * messages are written as loglevel increases. The 00076 * messages are written to a file in HTML format for viewing 00077 * in a web browser. @see HTML_logs 00078 * 00079 * @ingroup equilfunctions 00080 */ 00081 int vcs_equilibrate(thermo_t& s, const char* XY, 00082 int estimateEquil = 0, int printLvl = 0, 00083 int solver = -1, doublereal rtol = 1.0e-9, 00084 int maxsteps = 5000, 00085 int maxiter = 100, int loglevel = -99); 00086 00087 00088 //! Set a multi-phase chemical solution to chemical equilibrium. 00089 /*! 00090 * This function uses the vcs_MultiPhaseEquil interface to the 00091 * vcs solver. 00092 * The function uses the element abundance vector that is 00093 * currently consistent with the composition within the phases 00094 * themselves. Two other thermodynamic quantities, determined by the 00095 * XY string, are held constant during the equilibration. 00096 * 00097 * @param s The object to set to an equilibrium state 00098 * 00099 * @param XY A character string representing the unknowns 00100 * to be held constant 00101 * 00102 * @param estimateEquil integer indicating whether the solver 00103 * should estimate its own initial condition. 00104 * If 0, the initial mole fraction vector 00105 * in the %ThermoPhase object is used as the 00106 * initial condition. 00107 * If 1, the initial mole fraction vector 00108 * is used if the element abundances are 00109 * satisfied. 00110 * if -1, the initial mole fraction vector 00111 * is thrown out, and an estimate is 00112 * formulated. 00113 * 00114 * @param printLvl Determines the amount of printing that 00115 * gets sent to stdout from the vcs package 00116 * (Note, you may have to compile with debug 00117 * flags to get some printing). 00118 * 00119 * @param solver Determines which solver is used. 00120 * - 1 MultiPhaseEquil solver 00121 * - 2 VCSnonideal Solver (default) 00122 * 00123 * @param rtol Relative tolerance of the solve. Defaults to 00124 * 1.0E-9. 00125 * 00126 * @param maxsteps The maximum number of steps to take to find 00127 * the solution. 00128 * 00129 * @param maxiter For the MultiPhaseEquil solver only, this is 00130 * the maximum number of outer temperature or 00131 * pressure iterations to take when T and/or P is 00132 * not held fixed. 00133 * 00134 * @param loglevel Controls amount of diagnostic output. loglevel 00135 * = 0 suppresses diagnostics, and increasingly-verbose 00136 * messages are written as loglevel increases. The 00137 * messages are written to a file in HTML format for viewing 00138 * in a web browser. @see HTML_logs 00139 * 00140 * @ingroup equilfunctions 00141 */ 00142 int vcs_equilibrate(MultiPhase& s, const char* XY, 00143 int estimateEquil = 0, int printLvl = 0, 00144 int solver = 2, 00145 doublereal rtol = 1.0e-9, int maxsteps = 5000, 00146 int maxiter = 100, int loglevel = -99); 00147 00148 //! Set a multi-phase chemical solution to chemical equilibrium. 00149 /*! 00150 * This function uses the vcs_MultiPhaseEquil interface to the 00151 * vcs solver. 00152 * The function uses the element abundance vector that is 00153 * currently consistent with the composition within the phases 00154 * themselves. Two other thermodynamic quantities, determined by the 00155 * XY string, are held constant during the equilibration. 00156 * 00157 * @param s The object to set to an equilibrium state 00158 * 00159 * @param ixy An integer specifying the two properties to be held 00160 * constant. 00161 * 00162 * @param estimateEquil integer indicating whether the solver 00163 * should estimate its own initial condition. 00164 * If 0, the initial mole fraction vector 00165 * in the %ThermoPhase object is used as the 00166 * initial condition. 00167 * If 1, the initial mole fraction vector 00168 * is used if the element abundances are 00169 * satisfied. 00170 * if -1, the initial mole fraction vector 00171 * is thrown out, and an estimate is 00172 * formulated. 00173 * 00174 * @param printLvl Determines the amount of printing that 00175 * gets sent to stdout from the vcs package 00176 * (Note, you may have to compile with debug 00177 * flags to get some printing). 00178 * 00179 * @param solver Determines which solver is used. 00180 * - 1 MultiPhaseEquil solver 00181 * - 2 VCSnonideal Solver (default) 00182 * 00183 * @param rtol Relative tolerance of the solve. Defaults to 00184 * 1.0E-9. 00185 * 00186 * @param maxsteps The maximum number of steps to take to find 00187 * the solution. 00188 * 00189 * @param maxiter For the MultiPhaseEquil solver only, this is 00190 * the maximum number of outer temperature or 00191 * pressure iterations to take when T and/or P is 00192 * not held fixed. 00193 * 00194 * @param loglevel Controls amount of diagnostic output. loglevel 00195 * = 0 suppresses diagnostics, and increasingly-verbose 00196 * messages are written as loglevel increases. The 00197 * messages are written to a file in HTML format for viewing 00198 * in a web browser. @see HTML_logs 00199 * 00200 * @ingroup equilfunctions 00201 */ 00202 int vcs_equilibrate_1(MultiPhase& s, int ixy, 00203 int estimateEquil = 0, int printLvl = 0, 00204 int solver = 2, 00205 doublereal rtol = 1.0e-9, int maxsteps = 5000, 00206 int maxiter = 100, int loglevel = -99); 00207 00208 } 00209 00210 namespace VCSnonideal { 00211 00212 00213 class VCS_PROB; 00214 class VCS_SOLVE; 00215 00216 //! Translate a MultiPhase object into a VCS_PROB problem definition object 00217 /*! 00218 * @param mphase MultiPhase object that is the source for all of the information 00219 * @param vprob VCS_PROB problem definition that gets all of the information 00220 * 00221 * Note, both objects share the underlying Thermophase objects. So, neither 00222 * can be const objects. 00223 */ 00224 int vcs_Cantera_to_vprob(Cantera::MultiPhase *mphase, 00225 VCSnonideal::VCS_PROB *vprob); 00226 00227 //! Translate a MultiPhase information into a VCS_PROB problem definition object 00228 /*! 00229 * This version updates the problem statement information only. All species and 00230 * phase definitions remain the same. 00231 * 00232 * @param mphase MultiPhase object that is the source for all of the information 00233 * @param vprob VCS_PROB problem definition that gets all of the information 00234 * 00235 */ 00236 int vcs_Cantera_update_vprob(Cantera::MultiPhase *mphase, 00237 VCSnonideal::VCS_PROB *vprob); 00238 00239 //! Cantera's Interface to the Multiphase chemical equilibrium solver. 00240 /*! 00241 * Class MultiPhaseEquil is designed to be used to set a mixture 00242 * containing one or more phases to a state of chemical equilibrium. 00243 * 00244 * Note, as currently constructed, the underlying ThermoPhase 00245 * objects are shared between the MultiPhase object and this 00246 * object. Therefore, mix is not a const argument, and the 00247 * return parameters are contained in underlying ThermoPhase 00248 * objects. 00249 * 00250 * @ingroup equilfunctions 00251 */ 00252 class vcs_MultiPhaseEquil { 00253 public: 00254 00255 //! Shorthand for the MultiPhase mixture object used by Cantera 00256 //! to store information about multiple phases 00257 typedef Cantera::MultiPhase mix_t; 00258 00259 //! Typedef for an index variable 00260 typedef size_t index_t; 00261 00262 //! Typedef for a dense 2d matrix. 00263 typedef Cantera::DenseMatrix matrix_t; 00264 00265 //! Default empty constructor 00266 vcs_MultiPhaseEquil(); 00267 00268 00269 //! Constructor for the multiphase equilibrium solver 00270 /*! 00271 * This constructor will initialize the object with a MultiPhase 00272 * object, setting up the internal equilibration problem. 00273 * Note, as currently constructed, the underlying ThermoPhase 00274 * objects are shared between the MultiPhase object and this 00275 * object. Therefore, mix is not a const argument, and the 00276 * return parameters are contained in underlying ThermoPhase 00277 * objects. 00278 * 00279 * @param mix Object containing the MultiPhase object 00280 * @param printLvl Determines the amount of printing to stdout 00281 * that occurs for each call: 00282 * - 0 No printing 00283 * - 1 Only printing to the .csv file 00284 * - 2 print the soln only 00285 * - 3 Print the setup and then the soln only 00286 * - 4 Print a table for each iteration 00287 * - 5 Print more than a table for each iteration 00288 * 00289 */ 00290 vcs_MultiPhaseEquil(mix_t* mix, int printLvl); 00291 00292 //! Destructor for the class 00293 virtual ~vcs_MultiPhaseEquil(); 00294 00295 //! Return the index of the ith component 00296 /*! 00297 * Returns the index of the ith component in the equilibrium 00298 * calculation. The index refers to the ordering of the species 00299 * in the MultiPhase object. 00300 * 00301 * @param m Index of the component. Must be between 0 and the 00302 * number of components, which can be obtained from the 00303 * numComponents() command. 00304 */ 00305 int component(int m) const ; 00306 00307 //! Get the stoichiometric reaction coefficients for a single 00308 //! reaction index 00309 /*! 00310 * This returns a stoichiometric reaction vector for a single 00311 * formation reaction for a noncomponent species. There are 00312 * (nSpecies() - nComponents) formation reactions. Each 00313 * formation reaction will have a value of 1.0 for the species 00314 * that is being formed, and the other non-zero coefficients will 00315 * all involve the components of the mixture. 00316 * 00317 * @param rxn Reaction number. 00318 * @param nu Vector of coefficients for the formation reaction. 00319 * Length is equal to the number of species in 00320 * the MultiPhase object. 00321 */ 00322 void getStoichVector(index_t rxn, Cantera::vector_fp& nu); 00323 00324 //! return the number of iterations 00325 int iterations() const { return m_iter; } 00326 00327 //! Equilibrate the solution using the current element abundances 00328 //! storred in the MultiPhase object 00329 /*! 00330 * Use the vcs algorithm to equilibrate the current multiphase 00331 * mixture. 00332 * 00333 * @param XY Integer representing what two thermo quantities 00334 * are held constant during the equilibration 00335 * 00336 * @param estimateEquil integer indicating whether the solver 00337 * should estimate its own initial condition. 00338 * If 0, the initial mole fraction vector 00339 * in the %ThermoPhase object is used as the 00340 * initial condition. 00341 * If 1, the initial mole fraction vector 00342 * is used if the element abundances are 00343 * satisfied. 00344 * if -1, the initial mole fraction vector 00345 * is thrown out, and an estimate is 00346 * formulated. 00347 * 00348 * @param printLvl Determines the amount of printing that 00349 * gets sent to stdout from the vcs package 00350 * (Note, you may have to compile with debug 00351 * flags to get some printing). 00352 * @param err Internal error level 00353 * @param maxsteps max steps allowed. 00354 * @param loglevel for 00355 */ 00356 int equilibrate(int XY, int estimateEquil = 0, 00357 int printLvl= 0, doublereal err = 1.0e-6, 00358 int maxsteps = 5000, int loglevel=-99); 00359 00360 //! Equilibrate the solution using the current element abundances 00361 //! storred in the MultiPhase object using constant T and P 00362 /*! 00363 * Use the vcs algorithm to equilibrate the current multiphase 00364 * mixture. 00365 * 00366 * @param estimateEquil integer indicating whether the solver 00367 * should estimate its own initial condition. 00368 * If 0, the initial mole fraction vector 00369 * in the %ThermoPhase object is used as the 00370 * initial condition. 00371 * If 1, the initial mole fraction vector 00372 * is used if the element abundances are 00373 * satisfied. 00374 * if -1, the initial mole fraction vector 00375 * is thrown out, and an estimate is 00376 * formulated. 00377 * 00378 * @param printLvl Determines the amount of printing that 00379 * gets sent to stdout from the vcs package 00380 * (Note, you may have to compile with debug 00381 * flags to get some printing). 00382 * @param err Internal error level 00383 * @param maxsteps max steps allowed. 00384 * @param loglevel for 00385 */ 00386 int equilibrate_TP(int estimateEquil = 0, 00387 int printLvl= 0, doublereal err = 1.0e-6, 00388 int maxsteps = 5000, int loglevel=-99); 00389 00390 //! Equilibrate the solution using the current element abundances 00391 //! storred in the MultiPhase object using either constant H and P 00392 //! or constant U and P. 00393 /*! 00394 * Use the vcs algorithm to equilibrate the current multiphase 00395 * mixture. The pressure of the calculation is taken from 00396 * the current pressure storred with the MultiPhase object. 00397 * 00398 * @param Htarget Value of the total mixture enthalpy or total 00399 * internal energy that will be 00400 * kept constant. Note, this is and must be an extensive 00401 * quantity. units = Joules 00402 * 00403 * @param XY Integer flag indicating what is held constant. 00404 * Must be either HP or UP. 00405 * 00406 * @param Tlow Lower limit of the temperature. It's an 00407 * error condition if the temperature falls 00408 * below Tlow. 00409 * 00410 * @param Thigh Upper limit of the temperature. It's an 00411 * error condition if the temperature goes 00412 * higher than Thigh. 00413 * 00414 * @param estimateEquil integer indicating whether the solver 00415 * should estimate its own initial condition. 00416 * If 0, the initial mole fraction vector 00417 * in the %ThermoPhase object is used as the 00418 * initial condition. 00419 * If 1, the initial mole fraction vector 00420 * is used if the element abundances are 00421 * satisfied. 00422 * if -1, the initial mole fraction vector 00423 * is thrown out, and an estimate is 00424 * formulated. 00425 * 00426 * @param printLvl Determines the amount of printing that 00427 * gets sent to stdout from the vcs package 00428 * (Note, you may have to compile with debug 00429 * flags to get some printing). See main 00430 * constructor call for meaning of the levels. 00431 * 00432 * @param err Internal error level 00433 * 00434 * @param maxsteps max steps allowed. 00435 * 00436 * @param loglevel Determines the amount of printing to the HTML 00437 * output file. 00438 */ 00439 int equilibrate_HP(doublereal Htarget, int XY, double Tlow, double Thigh, 00440 int estimateEquil = 0, 00441 int printLvl = 0, doublereal err = 1.0E-6, 00442 int maxsteps = 5000, int loglevel=-99); 00443 00444 //! Equilibrate the solution using the current element abundances 00445 //! storred in the MultiPhase object using constant S and P. 00446 /*! 00447 * Use the vcs algorithm to equilibrate the current multiphase 00448 * mixture. The pressure of the calculation is taken from 00449 * the current pressure storred with the MultiPhase object. 00450 * 00451 * @param Starget Value of the total mixture entropy 00452 * that will be 00453 * kept constant. Note, this is and must be an extensive 00454 * quantity. units = Joules/K 00455 * 00456 * 00457 * @param Tlow Lower limit of the temperature. It's an 00458 * error condition if the temperature falls 00459 * below Tlow. 00460 * 00461 * @param Thigh Upper limit of the temperature. It's an 00462 * error condition if the temperature goes 00463 * higher than Thigh. 00464 * 00465 * @param estimateEquil integer indicating whether the solver 00466 * should estimate its own initial condition. 00467 * If 0, the initial mole fraction vector 00468 * in the %ThermoPhase object is used as the 00469 * initial condition. 00470 * If 1, the initial mole fraction vector 00471 * is used if the element abundances are 00472 * satisfied. 00473 * if -1, the initial mole fraction vector 00474 * is thrown out, and an estimate is 00475 * formulated. 00476 * 00477 * @param printLvl Determines the amount of printing that 00478 * gets sent to stdout from the vcs package 00479 * (Note, you may have to compile with debug 00480 * flags to get some printing). See main 00481 * constructor call for meaning of the levels. 00482 * 00483 * @param err Internal error level 00484 * 00485 * @param maxsteps max steps allowed. 00486 * 00487 * @param loglevel Determines the amount of printing to the HTML 00488 * output file. 00489 */ 00490 int equilibrate_SP(doublereal Starget, double Tlow, double Thigh, 00491 int estimateEquil = 0, 00492 int printLvl = 0, doublereal err = 1.0E-6, 00493 int maxsteps = 5000, int loglevel=-99); 00494 00495 00496 //! Equilibrate the solution using the current element abundances 00497 //! storred in the MultiPhase object using constant V and constant 00498 //! T, H, U, or S. 00499 /*! 00500 * Use the vcs algorithm to equilibrate the current multiphase 00501 * mixture. The pressure of the calculation is taken from 00502 * the current pressure storred with the MultiPhase object. 00503 * 00504 * 00505 * @param XY Integer flag indicating what is held constant. 00506 * Must be either TV, HV, UV, or SV. 00507 * 00508 * @param xtarget Value of the total thermodynamic parameter to 00509 * be held constant in addition to V. 00510 * Note, except for T, this must be an extensive 00511 * quantity. units = Joules/K or Joules 00512 * 00513 * @param estimateEquil integer indicating whether the solver 00514 * should estimate its own initial condition. 00515 * If 0, the initial mole fraction vector 00516 * in the %ThermoPhase object is used as the 00517 * initial condition. 00518 * If 1, the initial mole fraction vector 00519 * is used if the element abundances are 00520 * satisfied. 00521 * if -1, the initial mole fraction vector 00522 * is thrown out, and an estimate is 00523 * formulated. 00524 * 00525 * @param printLvl Determines the amount of printing that 00526 * gets sent to stdout from the vcs package 00527 * (Note, you may have to compile with debug 00528 * flags to get some printing). See main 00529 * constructor call for meaning of the levels. 00530 * 00531 * @param err Internal error level 00532 * 00533 * @param maxsteps max steps allowed. 00534 * 00535 * @param loglevel Determines the amount of printing to the HTML 00536 * output file. 00537 */ 00538 int equilibrate_TV(int XY, doublereal xtarget, 00539 int estimateEquil = 0, 00540 int printLvl = 0, doublereal err = 1.0E-6, 00541 int maxsteps = 5000, int loglevel = -99); 00542 00543 //! Report the equilibrium answer in a comma separated table format 00544 /*! 00545 * This routine is used for in the test suite. 00546 * 00547 * @param reportFile Base name of the file to get the report. 00548 * File name is incremented by 1 for each report. 00549 */ 00550 void reportCSV(const std::string &reportFile); 00551 00552 //! reports the number of components in the equilibration problem 00553 /*! 00554 * @return returns the number of components. If an equilibrium 00555 * problem hasn't been solved yet, it returns -1. 00556 */ 00557 int numComponents() const; 00558 00559 //! Reports the number of element contraints in the equilibration problem 00560 /*! 00561 * @return returns the number of element constraints. If an equilibrium 00562 * problem hasn't been solved yet, it returns -1. 00563 */ 00564 int numElemConstraints() const; 00565 00566 // Friend functions 00567 00568 friend int vcs_Cantera_to_vprob(Cantera::MultiPhase *mphase, 00569 VCSnonideal::VCS_PROB *vprob); 00570 friend int vcs_Cantera_update_vprob(Cantera::MultiPhase *mphase, 00571 VCSnonideal::VCS_PROB *vprob); 00572 00573 protected: 00574 00575 //! Vector that takes into account of the current sorting of the species 00576 /*! 00577 * The index of m_order is the original k value of the species in the 00578 * multiphase. The value of m_order, k_sorted, is the current value of the 00579 * species index. 00580 * 00581 * m_order[korig] = k_sorted 00582 */ 00583 Cantera::vector_int m_order; 00584 00585 //! Object which contains the problem statement 00586 /*! 00587 * The problem statement may contain some subtleties. For example, 00588 * the element constraints may be different than just an element 00589 * conservation contraint equations. 00590 * There may be kinetically frozen degrees of freedom. 00591 * There may be multiple electrolyte phases with zero charge constraints. 00592 * All of these make the problem statement different than the 00593 * simple element conservation statement. 00594 */ 00595 VCSnonideal::VCS_PROB *m_vprob; 00596 00597 //! Pointer to the MultiPhase mixture that will be equilibrated. 00598 /*! 00599 * Equilibrium solutions will be returned via this variable. 00600 */ 00601 mix_t *m_mix; 00602 00603 //! Print level from the VCSnonlinear package 00604 /*! 00605 * (Note, you may have to compile with debug 00606 * flags to get some printing). 00607 * 00608 * - 0 No IO from the routine whatsoever 00609 * - 1 file IO from reportCSV() carried out. 00610 * One line print statements from equilibrate_XY() functions 00611 * - 2 Problem statement information from vcs_Cantera_update_vprob() 00612 * - Final state of the system from vcs_solve_TP() 00613 * - 3 Several more setup tables 00614 * - Problem initialization routine 00615 * - 4 One table for each iteration within vcs_solve_Tp() 00616 * - 5 Multiple tables for each iteration within vcs_solve_TP() 00617 * - full discussion of decisions made for each variable. 00618 */ 00619 int m_printLvl; 00620 00621 //! Stoichiometric matrix 00622 /*! 00623 * 00624 */ 00625 matrix_t m_N; 00626 00627 //! Iteration Count 00628 int m_iter; 00629 00630 //! Vector of indices for species that are included in the 00631 //! calculation. 00632 /*! 00633 * This is used to exclude pure-phase species 00634 * with invalid thermo data 00635 */ 00636 Cantera::vector_int m_species; 00637 00638 //! Pointer to the object that does all of the equilibration work. 00639 /*! 00640 * VCS_SOLVE will have different ordering for species and element constraints 00641 * than this object or the VCS_PROB object. 00642 * This object owns the pointer. 00643 */ 00644 VCSnonideal::VCS_SOLVE *m_vsolvePtr; 00645 }; 00646 00647 } 00648 00649 00650 #endif
1.6.3