From 8522095deaacc3659b3195c5a5638bbe5b6aa606 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Wed, 16 Aug 2017 18:52:36 -0400 Subject: [PATCH] [Equil] Refactor to eliminate class VCS_PROB Move data in the VCS_PROB class to VCS_SOLVE --- include/cantera/equil/vcs_MultiPhaseEquil.h | 41 +-- include/cantera/equil/vcs_VolPhase.h | 3 +- include/cantera/equil/vcs_internal.h | 1 - include/cantera/equil/vcs_prob.h | 269 -------------- include/cantera/equil/vcs_solve.h | 253 +++++++++++-- src/equil/vcs_MultiPhaseEquil.cpp | 234 ++++++------ src/equil/vcs_prep.cpp | 8 +- src/equil/vcs_prob.cpp | 144 ++------ src/equil/vcs_solve.cpp | 384 +++++++++----------- 9 files changed, 550 insertions(+), 787 deletions(-) delete mode 100644 include/cantera/equil/vcs_prob.h diff --git a/include/cantera/equil/vcs_MultiPhaseEquil.h b/include/cantera/equil/vcs_MultiPhaseEquil.h index 2ca4e52a4..5f147f05f 100644 --- a/include/cantera/equil/vcs_MultiPhaseEquil.h +++ b/include/cantera/equil/vcs_MultiPhaseEquil.h @@ -11,30 +11,29 @@ #include "MultiPhase.h" #include "vcs_solve.h" -#include "vcs_prob.h" namespace Cantera { -//! Translate a MultiPhase object into a VCS_PROB problem definition object +//! Translate a MultiPhase object into a VCS_SOLVE problem definition object /*! * @param mphase MultiPhase object that is the source for all of the information - * @param vprob VCS_PROB problem definition that gets all of the information + * @param vsolve VCS_SOLVE problem definition that gets all of the information * * Note, both objects share the underlying ThermoPhase objects. So, neither can * be const objects. */ -int vcs_Cantera_to_vprob(MultiPhase* mphase, VCS_PROB* vprob); +int vcs_Cantera_to_vprob(MultiPhase* mphase, VCS_SOLVE* vsolve); -//! Translate a MultiPhase information into a VCS_PROB problem definition object +//! Translate a MultiPhase information into a VCS_SOLVE problem definition object /*! * This version updates the problem statement information only. All species and * phase definitions remain the same. * * @param mphase MultiPhase object that is the source for all of the information - * @param vprob VCS_PROB problem definition that gets all of the information + * @param vsolve VCS_SOLVE problem definition that gets all of the information */ -int vcs_Cantera_update_vprob(MultiPhase* mphase, VCS_PROB* vprob); +int vcs_Cantera_update_vprob(MultiPhase* mphase, VCS_SOLVE* vsolve); //! %Cantera's Interface to the Multiphase chemical equilibrium solver. /*! @@ -51,8 +50,6 @@ int vcs_Cantera_update_vprob(MultiPhase* mphase, VCS_PROB* vprob); class vcs_MultiPhaseEquil { public: - vcs_MultiPhaseEquil(); - //! Constructor for the multiphase equilibrium solver /*! * This constructor will initialize the object with a MultiPhase object, @@ -247,8 +244,8 @@ public: void reportCSV(const std::string& reportFile); // Friend functions - friend int vcs_Cantera_to_vprob(MultiPhase* mphase, VCS_PROB* vprob); - friend int vcs_Cantera_update_vprob(MultiPhase* mphase, VCS_PROB* vprob); + friend int vcs_Cantera_to_vprob(MultiPhase* mphase, VCS_SOLVE* vsolve); + friend int vcs_Cantera_update_vprob(MultiPhase* mphase, VCS_SOLVE* vsolve); protected: //! Vector that takes into account of the current sorting of the species @@ -261,17 +258,6 @@ protected: */ vector_int m_order; - //! Object which contains the problem statement - /*! - * The problem statement may contain some subtleties. For example, the - * element constraints may be different than just an element conservation - * contraint equations. There may be kinetically frozen degrees of freedom. - * There may be multiple electrolyte phases with zero charge constraints. - * All of these make the problem statement different than the simple element - * conservation statement. - */ - VCS_PROB m_vprob; - //! Pointer to the MultiPhase mixture that will be equilibrated. /*! * Equilibrium solutions will be returned via this variable. @@ -303,10 +289,17 @@ protected: //! is used to exclude pure-phase species with invalid thermo data vector_int m_species; - //! The object that does all of the equilibration work. + //! The object that contains the problem statement and does all of the equilibration work /*! + * The problem statement may contain some subtleties. For example, the + * element constraints may be different than just an element conservation + * contraint equations. There may be kinetically frozen degrees of freedom. + * There may be multiple electrolyte phases with zero charge constraints. + * All of these make the problem statement different than the simple element + * conservation statement. + * * VCS_SOLVE will have different ordering for species and element constraints - * than this object or the VCS_PROB object. + * than this object. */ VCS_SOLVE m_vsolve; }; diff --git a/include/cantera/equil/vcs_VolPhase.h b/include/cantera/equil/vcs_VolPhase.h index e70ecca48..2188a611c 100644 --- a/include/cantera/equil/vcs_VolPhase.h +++ b/include/cantera/equil/vcs_VolPhase.h @@ -622,8 +622,7 @@ private: */ vector_int m_speciesUnknownType; - //! Index of the element number in the global list of elements - //! stored in VCS_PROB or VCS_SOLVE + //! Index of the element number in the global list of elements stored in VCS_SOLVE std::vector m_elemGlobalIndex; //! Number of species in the phase diff --git a/include/cantera/equil/vcs_internal.h b/include/cantera/equil/vcs_internal.h index aefadead8..50b11d7f7 100644 --- a/include/cantera/equil/vcs_internal.h +++ b/include/cantera/equil/vcs_internal.h @@ -27,7 +27,6 @@ extern int vcs_timing_print_lvl; // Forward references class VCS_SPECIES_THERMO; -class VCS_PROB; //! Class to keep track of time and iterations /*! diff --git a/include/cantera/equil/vcs_prob.h b/include/cantera/equil/vcs_prob.h deleted file mode 100644 index 42bfb9fad..000000000 --- a/include/cantera/equil/vcs_prob.h +++ /dev/null @@ -1,269 +0,0 @@ -/** - * @file vcs_prob.h - * Header for the Interface class for the vcs thermo equilibrium solver package, - */ - -// This file is part of Cantera. See License.txt in the top-level directory or -// at http://www.cantera.org/license.txt for license and copyright information. - -#ifndef _VCS_PROB_H -#define _VCS_PROB_H - -#include "cantera/base/Array.h" - -namespace Cantera -{ - -class vcs_VolPhase; -class VCS_SPECIES_THERMO; - -//! Interface class for the vcs thermo equilibrium solver package, -//! which generally describes the problem to be solved. -class VCS_PROB -{ -public: - //! Problem type. I.e., the identity of what is held constant. Currently, T - //! and P are held constant, and this input is ignored - int prob_type; - - //! Total number of species in the problems - size_t nspecies; - - //! Species number used to size data structures - size_t NSPECIES0; - - //! Number of element constraints in the equilibrium problem - size_t ne; - - //! Number of element constraints used to size data structures - //! involving elements - size_t NE0; - - //! Number of phases in the problem - size_t NPhase; - - //! Number of phases used to size data structures - size_t NPHASE0; - - //! Vector of chemical potentials of the species. This is a calculated - //! output quantity. length = number of species. - vector_fp m_gibbsSpecies; - - //! Total number of moles of the kth species. - /*! - * This is both an input and an output variable. On input, this is an - * estimate of the mole numbers. The actual element abundance vector - * contains the problem specification. - * - * On output, this contains the solution for the total number of moles of - * the kth species. - */ - vector_fp w; - - //! Mole fraction vector. This is a calculated vector, calculated from w[]. - //! length number of species. - vector_fp mf; - - //! Element abundances for jth element - /*! - * This is input from the input file and is considered a constant from - * thereon within the vcs_solve_TP(). - */ - vector_fp gai; - - //! Formula Matrix for the problem - /*! - * FormulaMatrix(kspec,j) = Number of elements, j, in the kspec species - */ - Array2D FormulaMatrix; - - //! Specifies the species unknown type - /*! - * There are two types. One is the straightforward species, with the mole - * number w[k], as the unknown. The second is the an interfacial voltage - * where w[k] refers to the interfacial voltage in volts. - * - * These species types correspond to metallic electrons corresponding to - * electrodes. The voltage and other interfacial conditions sets up an - * interfacial current, which is set to zero in this initial treatment. - * Later we may have non-zero interfacial currents. - */ - vector_int SpeciesUnknownType; - - //! Temperature (Kelvin) - /*! - * Specification of the temperature for the equilibrium problem - */ - double T; - - //! Pressure - double PresPA; - - //! Volume of the entire system - /*! - * Note, this is an output variable atm - */ - double Vol; - - //! Partial Molar Volumes of species - /*! - * This is a calculated vector, calculated from w[]. - * length number of species. - */ - vector_fp VolPM; - - //! Specification of the initial estimate method - /*! - * * 0: user estimate - * * 1: user estimate if satisifies elements - * * -1: machine estimate - */ - int iest; - - //! Tolerance requirement for major species - double tolmaj; - - //! Tolerance requirement for minor species - double tolmin; - - //! Mapping between the species and the phases - std::vector PhaseID; - - //! Vector of strings containing the species names - std::vector SpName; - - //! vector of strings containing the element names - std::vector ElName; - - //! vector of Element types - vector_int m_elType; - - //! Specifies whether an element constraint is active - /*! - * The default is true - * Length = nelements - */ - vector_int ElActive; - - //! Molecular weight of species - /*! - * WtSpecies[k] = molecular weight of species in gm/mol - */ - vector_fp WtSpecies; - - //! Charge of each species - vector_fp Charge; - - //! Array of phase structures - std::vector VPhaseList; - - // String containing the title of the run - std::string Title; - - //! Vector of pointers to thermo structures which identify the model and - //! parameters for evaluating the thermodynamic functions for that - //! particular species - std::vector SpeciesThermo; - - //! Number of iterations. This is an output variable - int m_Iterations; - - //! Number of basis optimizations used. This is an output variable. - int m_NumBasisOptimizations; - - //! Print level for print routines - int m_printLvl; - - //! Debug print lvl - int vcs_debug_print_lvl; - - //! Constructor - /*! - * This constructor initializes the sizes within the object to parameter - * values. - * - * @param nsp number of species - * @param nel number of elements - * @param nph number of phases - */ - VCS_PROB(size_t nsp, size_t nel, size_t nph); - - ~VCS_PROB(); - - //! Resizes all of the element lists within the structure - /*! - * Note, this doesn't change the number of element constraints in the - * problem. It will change #NE0 if `nel` is greater than #NE0. - * - * @param nel size to dimension all the elements lists - * @param force If true, this will dimension the size to be equal to `nel` - * even if `nel` is less than the current value of #NE0 - */ - void resizeElements(size_t nel, int force); - - //! Calculate the element abundance vector from the mole numbers - void set_gai(); - - //! Print out the problem specification in all generality as it currently - //! exists in the VCS_PROB object - /*! - * @param print_lvl Parameter lvl for printing - * * 0 - no printing - * * 1 - all printing - */ - void prob_report(int print_lvl); - - //! Add elements to the local element list - /*! - * This routine sorts through the elements defined in the vcs_VolPhase - * object. It then adds the new elements to the VCS_PROB object, and creates - * a global map, which is stored in the vcs_VolPhase object. Id and matching - * of elements is done strictly via the element name, with case not - * mattering. - * - * The routine also fills in the position of the element in the vcs_VolPhase - * object's ElGlobalIndex field. - * - * @param volPhase Object containing the phase to be added. The elements in - * this phase are parsed for addition to the global element list - */ - void addPhaseElements(vcs_VolPhase* volPhase); - - //! This routine resizes the number of elements in the VCS_PROB object by - //! adding a new element to the end of the element list - /*! - * The element name is added. Formula vector entries ang element abundances - * for the new element are set to zero. - * - * @param elNameNew New name of the element - * @param elType Type of the element - * @param elactive boolean indicating whether the element is active - * @returns the index number of the new element - */ - size_t addElement(const char* elNameNew, int elType, int elactive); - - //! This routines adds entries for the formula matrix for one species - /*! - * This routines adds entries for the formula matrix for this object for one - * species - * - * This object also fills in the index filed, IndSpecies, within the - * volPhase object. - * - * @param volPhase object containing the species - * @param k Species number within the volPhase k - * @param kT global Species number within this object - * - */ - size_t addOnePhaseSpecies(vcs_VolPhase* volPhase, size_t k, size_t kT); - - //! Set the debug level - /*! - * @param vcs_debug_print_lvl input debug level - */ - void setDebugPrintLvl(int vcs_debug_print_lvl); -}; - -} - -#endif diff --git a/include/cantera/equil/vcs_solve.h b/include/cantera/equil/vcs_solve.h index 596ccc49b..0f1f13ac0 100644 --- a/include/cantera/equil/vcs_solve.h +++ b/include/cantera/equil/vcs_solve.h @@ -31,8 +31,8 @@ namespace Cantera class vcs_VolPhase; class VCS_SPECIES_THERMO; -class VCS_PROB; class VCS_COUNTERS; +class MultiPhase; //! This is the main structure used to hold the internal data //! used in vcs_solve_TP(), and to solve TP systems. @@ -44,10 +44,6 @@ class VCS_COUNTERS; class VCS_SOLVE { public: - VCS_SOLVE(); - - ~VCS_SOLVE(); - //! Initialize the sizes within the VCS_SOLVE object /*! * This resizes all of the internal arrays within the object. This routine @@ -63,20 +59,21 @@ public: * @param nelements Number of element constraints within the problem * @param nphase0 Number of phases defined within the problem. */ - void vcs_initSizes(const size_t nspecies0, const size_t nelements, const size_t nphase0); + VCS_SOLVE(size_t nspecies0, size_t nelements, size_t nphase0); + + ~VCS_SOLVE(); //! Solve an equilibrium problem /*! * This is the main interface routine to the equilibrium solver * - * @param vprob Object containing the equilibrium Problem statement * @param ifunc Determines the operation to be done: Valid values: * 0 -> Solve a new problem by initializing structures first. An * initial estimate may or may not have been already - * determined. This is indicated in the VCS_PROB structure. + * determined. This is indicated in the VCS_SOLVE structure. * 1 -> The problem has already been initialized and set up. We * call this routine to resolve it using the problem - * statement and solution estimate contained in the VCS_PROB + * statement and solution estimate contained in the VCS_SOLVE * structure. * 2 -> Don't solve a problem. Destroy all the private * structures. @@ -90,7 +87,7 @@ public: * @return nonzero value: failure to solve the problem at hand. zero : * success */ - int vcs(VCS_PROB* vprob, int ifunc, int ipr, int ip1, int maxit); + int vcs(int ifunc, int ipr, int ip1, int maxit); //! Main routine that solves for equilibrium at constant T and P using a //! variant of the VCS method @@ -587,12 +584,10 @@ public: * Current there is one condition. If all the element abundances are zero, * the algorithm will fail. * - * @param vprob VCS_PROB pointer to the definition of the equilibrium - * problem * @return If true, the problem is well-posed. If false, the problem * is not well posed. */ - bool vcs_wellPosed(VCS_PROB* vprob); + bool vcs_wellPosed(); //! Rearrange the constraint equations represented by the Formula //! Matrix so that the operational ones are in the front @@ -816,36 +811,17 @@ public: double vcs_GibbsPhase(size_t iphase, const double* const w, const double* const fe); - //! Transfer the results of the equilibrium calculation back to VCS_PROB - /*! - * The VCS_PROB structure is returned to the user. - * - * @param pub Pointer to VCS_PROB object that will get the results of the - * equilibrium calculation transfered to it. - */ - int vcs_prob_update(VCS_PROB* pub); + //! Transfer the results of the equilibrium calculation back from VCS_SOLVE + int vcs_prob_update(); - //! Fully specify the problem to be solved using VCS_PROB - /*! - * Use the contents of the VCS_PROB to specify the contents of the - * private data, VCS_SOLVE. - * - * @param pub Pointer to VCS_PROB that will be used to - * initialize the current equilibrium problem - */ - int vcs_prob_specifyFully(const VCS_PROB* pub); + //! Fully specify the problem to be solved + int vcs_prob_specifyFully(); - //! Specify the problem to be solved using VCS_PROB, incrementally + //! Specify the problem to be solved, incrementally /*! - * Use the contents of the VCS_PROB to specify the contents of the - * private data, VCS_SOLVE. - * * It's assumed we are solving the same problem. - * - * @param pub Pointer to VCS_PROB that will be used to initialize the - * current equilibrium problem */ - int vcs_prob_specify(const VCS_PROB* pub); + int vcs_prob_specify(); private: //! Zero out the concentration of a species. @@ -1085,6 +1061,203 @@ private: vector_fp m_wx; public: + //! @{ Variables moved from VCS_PROB + + //! Problem type. I.e., the identity of what is held constant. Currently, T + //! and P are held constant, and this input is ignored + int prob_type; + + //! Total number of species in the problems + size_t nspecies; + + //! Number of element constraints in the equilibrium problem + size_t ne; + //! Number of element constraints used to size data structures + //! involving elements + size_t NE0; + //! Number of phases in the problem + size_t NPhase; + //! Vector of chemical potentials of the species. This is a calculated + //! output quantity. length = number of species. + vector_fp m_gibbsSpecies; + //! Total number of moles of the kth species. + /*! + * This is both an input and an output variable. On input, this is an + * estimate of the mole numbers. The actual element abundance vector + * contains the problem specification. + * + * On output, this contains the solution for the total number of moles of + * the kth species. + */ + vector_fp w; + //! Mole fraction vector. This is a calculated vector, calculated from w[]. + //! length number of species. + vector_fp mf; + //! Element abundances for jth element + /*! + * This is input from the input file and is considered a constant from + * thereon within the vcs_solve_TP(). + */ + vector_fp gai; + + //! Formula Matrix for the problem + /*! + * FormulaMatrix(kspec,j) = Number of elements, j, in the kspec species + */ + Array2D FormulaMatrix; + + //! Specifies the species unknown type + /*! + * There are two types. One is the straightforward species, with the mole + * number w[k], as the unknown. The second is the an interfacial voltage + * where w[k] refers to the interfacial voltage in volts. + * + * These species types correspond to metallic electrons corresponding to + * electrodes. The voltage and other interfacial conditions sets up an + * interfacial current, which is set to zero in this initial treatment. + * Later we may have non-zero interfacial currents. + */ + vector_int SpeciesUnknownType; + + //! Temperature (Kelvin) + /*! + * Specification of the temperature for the equilibrium problem + */ + double T; + + //! Pressure + double PresPA; + + //! Volume of the entire system + /*! + * Note, this is an output variable atm + */ + double Vol; + + //! Partial Molar Volumes of species + /*! + * This is a calculated vector, calculated from w[]. + * length number of species. + */ + vector_fp VolPM; + + //! Specification of the initial estimate method + /*! + * * 0: user estimate + * * 1: user estimate if satisifies elements + * * -1: machine estimate + */ + int iest; + + //! Tolerance requirement for major species + double tolmaj; + + //! Tolerance requirement for minor species + double tolmin; + + //! Mapping between the species and the phases + std::vector PhaseID; + + //! Specifies whether an element constraint is active + /*! + * The default is true + * Length = nelements + */ + vector_int ElActive; + + //! Molecular weight of species + /*! + * WtSpecies[k] = molecular weight of species in gm/mol + */ + vector_fp WtSpecies; + + //! Charge of each species + vector_fp Charge; + + //! Array of phase structures + std::vector VPhaseList; + + // String containing the title of the run + std::string Title; + + //! Vector of pointers to thermo structures which identify the model and + //! parameters for evaluating the thermodynamic functions for that + //! particular species + std::vector SpeciesThermo; + + //! Number of iterations. This is an output variable + int m_Iterations; + + //! Number of basis optimizations used. This is an output variable. + int m_NumBasisOptimizations; + + //! Print level for print routines + int m_printLvl; + + //! Debug print lvl + int vcs_debug_print_lvl; + + MultiPhase* m_mix; + + //! Print out the problem specification in all generality as it currently + //! exists in the VCS_SOLVE object + /*! + * @param print_lvl Parameter lvl for printing + * * 0 - no printing + * * 1 - all printing + */ + void prob_report(int print_lvl); + + //! Add elements to the local element list + /*! + * This routine sorts through the elements defined in the vcs_VolPhase + * object. It then adds the new elements to the VCS_SOLVE object, and creates + * a global map, which is stored in the vcs_VolPhase object. Id and matching + * of elements is done strictly via the element name, with case not + * mattering. + * + * The routine also fills in the position of the element in the vcs_VolPhase + * object's ElGlobalIndex field. + * + * @param volPhase Object containing the phase to be added. The elements in + * this phase are parsed for addition to the global element list + */ + void addPhaseElements(vcs_VolPhase* volPhase); + + //! This routines adds entries for the formula matrix for one species + /*! + * This routines adds entries for the formula matrix for this object for one + * species + * + * This object also fills in the index filed, IndSpecies, within the + * volPhase object. + * + * @param volPhase object containing the species + * @param k Species number within the volPhase k + * @param kT global Species number within this object + * + */ + size_t addOnePhaseSpecies(vcs_VolPhase* volPhase, size_t k, size_t kT); + //! @} + + //! Calculate the element abundance vector from the mole numbers + void set_gai(); + + //! This routine resizes the number of elements in the VCS_SOLVE object by + //! adding a new element to the end of the element list + /*! + * The element name is added. Formula vector entries ang element abundances + * for the new element are set to zero. + * + * @param elNameNew New name of the element + * @param elType Type of the element + * @param elactive boolean indicating whether the element is active + * @returns the index number of the new element + */ + size_t addElement(const char* elNameNew, int elType, int elactive); + + void reportCSV(const std::string& reportFile); + //! value of the number of species used to size data structures size_t NSPECIES0; @@ -1355,7 +1528,7 @@ public: * k = m_speciesMapIndex[kspec] * * kspec = current order in the vcs_solve object - * k = original order in the vcs_prob object and in the MultiPhase object + * k = original order in the MultiPhase object */ std::vector m_speciesMapIndex; @@ -1382,7 +1555,7 @@ public: * * e = m_elementMapIndex[eNum] * eNum = current order in the vcs_solve object - * e = original order in the vcs_prob object and in the MultiPhase object + * e = original order in the MultiPhase object */ std::vector m_elementMapIndex; diff --git a/src/equil/vcs_MultiPhaseEquil.cpp b/src/equil/vcs_MultiPhaseEquil.cpp index 54e119b08..dbadb8887 100644 --- a/src/equil/vcs_MultiPhaseEquil.cpp +++ b/src/equil/vcs_MultiPhaseEquil.cpp @@ -21,24 +21,18 @@ using namespace std; namespace Cantera { -vcs_MultiPhaseEquil::vcs_MultiPhaseEquil() : - m_vprob(0, 0, 0), - m_mix(0), - m_printLvl(0) -{ -} - vcs_MultiPhaseEquil::vcs_MultiPhaseEquil(MultiPhase* mix, int printLvl) : - m_vprob(mix->nSpecies(), mix->nElements(), mix->nPhases()), m_mix(0), - m_printLvl(printLvl) + m_printLvl(printLvl), + m_vsolve(mix->nSpecies(), mix->nElements(), mix->nPhases()) { m_mix = mix; - m_vprob.m_printLvl = m_printLvl; + m_vsolve.m_printLvl = m_printLvl; + m_vsolve.m_mix = m_mix; - // Work out the details of the VCS_VPROB construction and Transfer the - // current problem to VCS_PROB object - int res = vcs_Cantera_to_vprob(mix, &m_vprob); + // Work out the details of the VCS_SOLVE construction and Transfer the + // current problem to the VCS_SOLVE object + int res = vcs_Cantera_to_vprob(mix, &m_vsolve); if (res != 0) { plogf("problems\n"); } @@ -435,20 +429,20 @@ int vcs_MultiPhaseEquil::equilibrate_TP(int estimateEquil, int maxit = maxsteps; clockWC tickTock; m_printLvl = printLvl; - m_vprob.m_printLvl = printLvl; + m_vsolve.m_printLvl = printLvl; // Extract the current state information from the MultiPhase object and // Transfer it to VCS_PROB object. - int res = vcs_Cantera_update_vprob(m_mix, &m_vprob); + int res = vcs_Cantera_update_vprob(m_mix, &m_vsolve); if (res != 0) { plogf("problems\n"); } // Set the estimation technique if (estimateEquil) { - m_vprob.iest = estimateEquil; + m_vsolve.iest = estimateEquil; } else { - m_vprob.iest = 0; + m_vsolve.iest = 0; } // Check obvious bounds on the temperature and pressure NOTE, we may want to @@ -464,7 +458,7 @@ int vcs_MultiPhaseEquil::equilibrate_TP(int estimateEquil, // Print out the problem specification from the point of // view of the vprob object. - m_vprob.prob_report(m_printLvl); + m_vsolve.prob_report(m_printLvl); //! Call the thermo Program int ip1 = m_printLvl; @@ -474,7 +468,7 @@ int vcs_MultiPhaseEquil::equilibrate_TP(int estimateEquil, } else { ip1 = 0; } - int iSuccess = m_vsolve.vcs(&m_vprob, 0, ipr, ip1, maxit); + int iSuccess = m_vsolve.vcs(0, ipr, ip1, maxit); // Transfer the information back to the MultiPhase object. Note we don't // just call setMoles, because some multispecies solution phases may be @@ -483,11 +477,11 @@ int vcs_MultiPhaseEquil::equilibrate_TP(int estimateEquil, // about likely reemergent states. m_mix->uploadMoleFractionsFromPhases(); size_t kGlob = 0; - for (size_t ip = 0; ip < m_vprob.NPhase; ip++) { + for (size_t ip = 0; ip < m_vsolve.NPhase; ip++) { double phaseMole = 0.0; ThermoPhase& tref = m_mix->phase(ip); for (size_t k = 0; k < tref.nSpecies(); k++, kGlob++) { - phaseMole += m_vprob.w[kGlob]; + phaseMole += m_vsolve.w[kGlob]; } m_mix->setPhaseMoles(ip, phaseMole); } @@ -499,8 +493,8 @@ int vcs_MultiPhaseEquil::equilibrate_TP(int estimateEquil, plogf("\nVCS FAILED TO CONVERGE!\n"); } plogf("\n"); - plogf("Temperature = %g Kelvin\n", m_vprob.T); - plogf("Pressure = %g Pa\n", m_vprob.PresPA); + plogf("Temperature = %g Kelvin\n", m_vsolve.T); + plogf("Pressure = %g Pa\n", m_vsolve.PresPA); plogf("\n"); plogf("----------------------------------------" "---------------------\n"); @@ -508,23 +502,23 @@ int vcs_MultiPhaseEquil::equilibrate_TP(int estimateEquil, plogf(" Mole_Fraction Chem_Potential (J/kmol)\n"); plogf("--------------------------------------------------" "-----------\n"); - for (size_t i = 0; i < m_vprob.nspecies; i++) { - plogf("%-12s", m_vprob.SpName[i]); - if (m_vprob.SpeciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - plogf(" %15.3e %15.3e ", 0.0, m_vprob.mf[i]); - plogf("%15.3e\n", m_vprob.m_gibbsSpecies[i]); + for (size_t i = 0; i < m_vsolve.nspecies; i++) { + plogf("%-12s", m_mix->speciesName(i)); + if (m_vsolve.SpeciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { + plogf(" %15.3e %15.3e ", 0.0, m_vsolve.mf[i]); + plogf("%15.3e\n", m_vsolve.m_gibbsSpecies[i]); } else { - plogf(" %15.3e %15.3e ", m_vprob.w[i], m_vprob.mf[i]); - if (m_vprob.w[i] <= 0.0) { - size_t iph = m_vprob.PhaseID[i]; - vcs_VolPhase* VPhase = m_vprob.VPhaseList[iph]; + plogf(" %15.3e %15.3e ", m_vsolve.w[i], m_vsolve.mf[i]); + if (m_vsolve.w[i] <= 0.0) { + size_t iph = m_vsolve.PhaseID[i]; + vcs_VolPhase* VPhase = m_vsolve.VPhaseList[iph]; if (VPhase->nSpecies() > 1) { plogf(" -1.000e+300\n"); } else { - plogf("%15.3e\n", m_vprob.m_gibbsSpecies[i]); + plogf("%15.3e\n", m_vsolve.m_gibbsSpecies[i]); } } else { - plogf("%15.3e\n", m_vprob.m_gibbsSpecies[i]); + plogf("%15.3e\n", m_vsolve.m_gibbsSpecies[i]); } } } @@ -539,15 +533,15 @@ int vcs_MultiPhaseEquil::equilibrate_TP(int estimateEquil, void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile) { - size_t nphase = m_vprob.NPhase; + size_t nphase = m_vsolve.NPhase; FILE* FP = fopen(reportFile.c_str(), "w"); if (!FP) { throw CanteraError("vcs_MultiPhaseEquil::reportCSV", "Failure to open file"); } - vector_fp& mf = m_vprob.mf; - double* fe = &m_vprob.m_gibbsSpecies[0]; + vector_fp& mf = m_vsolve.mf; + double* fe = &m_vsolve.m_gibbsSpecies[0]; vector_fp VolPM; vector_fp activity; vector_fp ac; @@ -562,7 +556,7 @@ void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile) size_t nSpecies = tref.nSpecies(); VolPM.resize(nSpecies, 0.0); tref.getPartialMolarVolumes(&VolPM[0]); - vcs_VolPhase* volP = m_vprob.VPhaseList[iphase]; + vcs_VolPhase* volP = m_vsolve.VPhaseList[iphase]; double TMolesPhase = volP->totalMoles(); double VolPhaseVolumes = 0.0; @@ -578,14 +572,14 @@ void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile) fprintf(FP,"Temperature = %11.5g kelvin\n", m_mix->temperature()); fprintf(FP,"Pressure = %11.5g Pascal\n", m_mix->pressure()); fprintf(FP,"Total Volume = %11.5g m**3\n", vol); - fprintf(FP,"Number Basis optimizations = %d\n", m_vprob.m_NumBasisOptimizations); - fprintf(FP,"Number VCS iterations = %d\n", m_vprob.m_Iterations); + fprintf(FP,"Number Basis optimizations = %d\n", m_vsolve.m_NumBasisOptimizations); + fprintf(FP,"Number VCS iterations = %d\n", m_vsolve.m_Iterations); for (size_t iphase = 0; iphase < nphase; iphase++) { size_t istart = m_mix->speciesIndex(0, iphase); ThermoPhase& tref = m_mix->phase(iphase); string phaseName = tref.name(); - vcs_VolPhase* volP = m_vprob.VPhaseList[iphase]; + vcs_VolPhase* volP = m_vsolve.VPhaseList[iphase]; double TMolesPhase = volP->totalMoles(); size_t nSpecies = tref.nSpecies(); activity.resize(nSpecies, 0.0); @@ -674,7 +668,7 @@ void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile) // HKM -> Work on transferring the current value of the voltages into the // equilibrium problem. -int vcs_Cantera_to_vprob(MultiPhase* mphase, VCS_PROB* vprob) +int vcs_Cantera_to_vprob(MultiPhase* mphase, VCS_SOLVE* vsolve) { VCS_SPECIES_THERMO* ts_ptr = 0; @@ -683,19 +677,19 @@ int vcs_Cantera_to_vprob(MultiPhase* mphase, VCS_PROB* vprob) size_t totNumSpecies = mphase->nSpecies(); // Problem type has yet to be worked out. - vprob->prob_type = 0; - vprob->nspecies = totNumSpecies; - vprob->ne = 0; - vprob->NPhase = totNumPhases; + vsolve->prob_type = 0; + vsolve->nspecies = totNumSpecies; + vsolve->ne = 0; + vsolve->NPhase = totNumPhases; // Set the initial estimate to a machine generated estimate for now // We will work out the details later. - vprob->iest = -1; - vprob->T = mphase->temperature(); - vprob->PresPA = mphase->pressure(); - vprob->Vol = mphase->volume(); - vprob->Title = "MultiPhase Object"; + vsolve->iest = -1; + vsolve->T = mphase->temperature(); + vsolve->PresPA = mphase->pressure(); + vsolve->Vol = mphase->volume(); + vsolve->Title = "MultiPhase Object"; - int printLvl = vprob->m_printLvl; + int printLvl = vsolve->m_printLvl; // Loop over the phases, transferring pertinent information int kT = 0; @@ -720,7 +714,7 @@ int vcs_Cantera_to_vprob(MultiPhase* mphase, VCS_PROB* vprob) // ->NumSpecies = number of species in the phase // ->TMolesInert = Inerts in the phase = 0.0 for cantera // ->PhaseName = Name of the phase - vcs_VolPhase* VolPhase = vprob->VPhaseList[iphase]; + vcs_VolPhase* VolPhase = vsolve->VPhaseList[iphase]; VolPhase->resize(iphase, nSpPhase, nelem, phaseName.c_str(), 0.0); VolPhase->m_gasPhase = gasPhase; @@ -765,8 +759,8 @@ int vcs_Cantera_to_vprob(MultiPhase* mphase, VCS_PROB* vprob) // Combine the element information in the vcs_VolPhase // object into the vprob object. - vprob->addPhaseElements(VolPhase); - VolPhase->setState_TP(vprob->T, vprob->PresPA); + vsolve->addPhaseElements(VolPhase); + VolPhase->setState_TP(vsolve->T, vsolve->PresPA); vector_fp muPhase(tPhase->nSpecies(),0.0); tPhase->getChemPotentials(&muPhase[0]); double tMoles = 0.0; @@ -775,63 +769,59 @@ int vcs_Cantera_to_vprob(MultiPhase* mphase, VCS_PROB* vprob) for (size_t k = 0; k < nSpPhase; k++) { // Obtain the molecular weight of the species from the // ThermoPhase object - vprob->WtSpecies[kT] = tPhase->molecularWeight(k); + vsolve->WtSpecies[kT] = tPhase->molecularWeight(k); // Obtain the charges of the species from the ThermoPhase object - vprob->Charge[kT] = tPhase->charge(k); + vsolve->Charge[kT] = tPhase->charge(k); // Set the phaseid of the species - vprob->PhaseID[kT] = iphase; - - // Transfer the Species name - string stmp = mphase->speciesName(kT); - vprob->SpName[kT] = stmp; + vsolve->PhaseID[kT] = iphase; // Transfer the type of unknown - vprob->SpeciesUnknownType[kT] = VolPhase->speciesUnknownType(k); - if (vprob->SpeciesUnknownType[kT] == VCS_SPECIES_TYPE_MOLNUM) { + vsolve->SpeciesUnknownType[kT] = VolPhase->speciesUnknownType(k); + if (vsolve->SpeciesUnknownType[kT] == VCS_SPECIES_TYPE_MOLNUM) { // Set the initial number of kmoles of the species // and the mole fraction vector - vprob->w[kT] = mphase->speciesMoles(kT); - tMoles += vprob->w[kT]; - vprob->mf[kT] = mphase->moleFraction(kT); - } else if (vprob->SpeciesUnknownType[kT] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - vprob->w[kT] = tPhase->electricPotential(); - vprob->mf[kT] = mphase->moleFraction(kT); + vsolve->w[kT] = mphase->speciesMoles(kT); + tMoles += vsolve->w[kT]; + vsolve->mf[kT] = mphase->moleFraction(kT); + } else if (vsolve->SpeciesUnknownType[kT] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { + vsolve->w[kT] = tPhase->electricPotential(); + vsolve->mf[kT] = mphase->moleFraction(kT); } else { - throw CanteraError(" vcs_Cantera_to_vprob() ERROR", - "Unknown species type: {}", vprob->SpeciesUnknownType[kT]); + throw CanteraError(" vcs_Cantera_to_vsolve() ERROR", + "Unknown species type: {}", vsolve->SpeciesUnknownType[kT]); } // transfer chemical potential vector - vprob->m_gibbsSpecies[kT] = muPhase[k]; + vsolve->m_gibbsSpecies[kT] = muPhase[k]; // Transfer the species information from the // volPhase structure to the VPROB structure // This includes: // FormulaMatrix[][] // VolPhase->IndSpecies[] - vprob->addOnePhaseSpecies(VolPhase, k, kT); + vsolve->addOnePhaseSpecies(VolPhase, k, kT); // Get a pointer to the thermo object - ts_ptr = vprob->SpeciesThermo[kT]; + ts_ptr = vsolve->SpeciesThermo[kT]; // Fill in the vcs_SpeciesProperty structure vcs_SpeciesProperties* sProp = VolPhase->speciesProperty(k); - sProp->NumElements = vprob->ne; - sProp->SpName = vprob->SpName[kT]; + sProp->NumElements = vsolve->ne; + sProp->SpName = mphase->speciesName(kT); sProp->SpeciesThermo = ts_ptr; sProp->WtSpecies = tPhase->molecularWeight(k); - sProp->FormulaMatrixCol.resize(vprob->ne, 0.0); - for (size_t e = 0; e < vprob->ne; e++) { - sProp->FormulaMatrixCol[e] = vprob->FormulaMatrix(kT,e); + sProp->FormulaMatrixCol.resize(vsolve->ne, 0.0); + for (size_t e = 0; e < vsolve->ne; e++) { + sProp->FormulaMatrixCol[e] = vsolve->FormulaMatrix(kT,e); } sProp->Charge = tPhase->charge(k); sProp->SurfaceSpecies = false; sProp->VolPM = 0.0; // Transfer the thermo specification of the species - // vprob->SpeciesThermo[] + // vsolve->SpeciesThermo[] // Add lookback connectivity into the thermo object first ts_ptr->IndexPhase = iphase; @@ -859,7 +849,7 @@ int vcs_Cantera_to_vprob(MultiPhase* mphase, VCS_PROB* vprob) ts_ptr->SSStar_Vol_Model = VCS_SSVOL_CONSTANT; } } else { - if (vprob->m_printLvl > 2) { + if (vsolve->m_printLvl > 2) { plogf("vcs_Cantera_convert: Species Type %d not known \n", spType); } @@ -887,18 +877,18 @@ int vcs_Cantera_to_vprob(MultiPhase* mphase, VCS_PROB* vprob) if (tMoles > 0.0) { for (size_t k = 0; k < nSpPhase; k++) { size_t kTa = VolPhase->spGlobalIndexVCS(k); - vprob->mf[kTa] = vprob->w[kTa] / tMoles; + vsolve->mf[kTa] = vsolve->w[kTa] / tMoles; } } else { // Perhaps, we could do a more sophisticated treatment below. // But, will start with this. for (size_t k = 0; k < nSpPhase; k++) { size_t kTa = VolPhase->spGlobalIndexVCS(k); - vprob->mf[kTa]= 1.0 / (double) nSpPhase; + vsolve->mf[kTa]= 1.0 / (double) nSpPhase; } } - VolPhase->setMolesFromVCS(VCS_STATECALC_OLD, &vprob->w[0]); + VolPhase->setMolesFromVCS(VCS_STATECALC_OLD, &vsolve->w[0]); // Now, calculate a sample naught Gibbs free energy calculation // at the specified temperature. @@ -906,17 +896,17 @@ int vcs_Cantera_to_vprob(MultiPhase* mphase, VCS_PROB* vprob) vcs_SpeciesProperties* sProp = VolPhase->speciesProperty(k); ts_ptr = sProp->SpeciesThermo; ts_ptr->SS0_feSave = VolPhase->G0_calc_one(k)/ GasConstant; - ts_ptr->SS0_TSave = vprob->T; + ts_ptr->SS0_TSave = vsolve->T; } } // Transfer initial element abundances to the vprob object. // We have to find the mapping index from one to the other - vprob->gai.resize(vprob->ne, 0.0); - vprob->set_gai(); + vsolve->gai.resize(vsolve->ne, 0.0); + vsolve->set_gai(); // Printout the species information: PhaseID's and mole nums - if (vprob->m_printLvl > 1) { + if (vsolve->m_printLvl > 1) { writeline('=', 80, true, true); writeline('=', 16, false); plogf(" Cantera_to_vprob: START OF PROBLEM STATEMENT "); @@ -925,16 +915,16 @@ int vcs_Cantera_to_vprob(MultiPhase* mphase, VCS_PROB* vprob) plogf(" Phase IDs of species\n"); plogf(" species phaseID phaseName "); plogf(" Initial_Estimated_kMols\n"); - for (size_t i = 0; i < vprob->nspecies; i++) { - size_t iphase = vprob->PhaseID[i]; + for (size_t i = 0; i < vsolve->nspecies; i++) { + size_t iphase = vsolve->PhaseID[i]; - vcs_VolPhase* VolPhase = vprob->VPhaseList[iphase]; - plogf("%16s %5d %16s", vprob->SpName[i].c_str(), iphase, + vcs_VolPhase* VolPhase = vsolve->VPhaseList[iphase]; + plogf("%16s %5d %16s", mphase->speciesName(i).c_str(), iphase, VolPhase->PhaseName.c_str()); - if (vprob->SpeciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - plogf(" Volts = %-10.5g\n", vprob->w[i]); + if (vsolve->SpeciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { + plogf(" Volts = %-10.5g\n", vsolve->w[i]); } else { - plogf(" %-10.5g\n", vprob->w[i]); + plogf(" %-10.5g\n", vsolve->w[i]); } } @@ -944,8 +934,8 @@ int vcs_Cantera_to_vprob(MultiPhase* mphase, VCS_PROB* vprob) plogf(" PhaseName PhaseNum SingSpec GasPhase EqnState NumSpec"); plogf(" TMolesInert Tmoles(kmol)\n"); - for (size_t iphase = 0; iphase < vprob->NPhase; iphase++) { - vcs_VolPhase* VolPhase = vprob->VPhaseList[iphase]; + for (size_t iphase = 0; iphase < vsolve->NPhase; iphase++) { + vcs_VolPhase* VolPhase = vsolve->VPhaseList[iphase]; plogf("%16s %5d %5d %8d %16s %8d %16e ", VolPhase->PhaseName.c_str(), VolPhase->VP_ID_, VolPhase->m_singleSpecies, VolPhase->m_gasPhase, VolPhase->eos_name(), @@ -963,29 +953,29 @@ int vcs_Cantera_to_vprob(MultiPhase* mphase, VCS_PROB* vprob) return VCS_SUCCESS; } -int vcs_Cantera_update_vprob(MultiPhase* mphase, VCS_PROB* vprob) +int vcs_Cantera_update_vprob(MultiPhase* mphase, VCS_SOLVE* vsolve) { size_t totNumPhases = mphase->nPhases(); size_t kT = 0; vector_fp tmpMoles; // Problem type has yet to be worked out. - vprob->prob_type = 0; + vsolve->prob_type = 0; // Whether we have an estimate or not gets overwritten on // the call to the equilibrium solver. - vprob->iest = -1; - vprob->T = mphase->temperature(); - vprob->PresPA = mphase->pressure(); - vprob->Vol = mphase->volume(); + vsolve->iest = -1; + vsolve->T = mphase->temperature(); + vsolve->PresPA = mphase->pressure(); + vsolve->Vol = mphase->volume(); for (size_t iphase = 0; iphase < totNumPhases; iphase++) { ThermoPhase* tPhase = &mphase->phase(iphase); - vcs_VolPhase* volPhase = vprob->VPhaseList[iphase]; + vcs_VolPhase* volPhase = vsolve->VPhaseList[iphase]; // Set the electric potential of the volume phase from the // ThermoPhase object's value. volPhase->setElectricPotential(tPhase->electricPotential()); - volPhase->setState_TP(vprob->T, vprob->PresPA); + volPhase->setState_TP(vsolve->T, vsolve->PresPA); vector_fp muPhase(tPhase->nSpecies(),0.0); tPhase->getChemPotentials(&muPhase[0]); @@ -994,20 +984,20 @@ int vcs_Cantera_update_vprob(MultiPhase* mphase, VCS_PROB* vprob) tmpMoles.resize(nSpPhase); for (size_t k = 0; k < nSpPhase; k++) { tmpMoles[k] = mphase->speciesMoles(kT); - vprob->w[kT] = mphase->speciesMoles(kT); - vprob->mf[kT] = mphase->moleFraction(kT); + vsolve->w[kT] = mphase->speciesMoles(kT); + vsolve->mf[kT] = mphase->moleFraction(kT); // transfer chemical potential vector - vprob->m_gibbsSpecies[kT] = muPhase[k]; + vsolve->m_gibbsSpecies[kT] = muPhase[k]; kT++; } if (volPhase->phiVarIndex() != npos) { size_t kphi = volPhase->phiVarIndex(); size_t kglob = volPhase->spGlobalIndexVCS(kphi); - vprob->w[kglob] = tPhase->electricPotential(); + vsolve->w[kglob] = tPhase->electricPotential(); } - volPhase->setMolesFromVCS(VCS_STATECALC_OLD, &vprob->w[0]); + volPhase->setMolesFromVCS(VCS_STATECALC_OLD, &vsolve->w[0]); if ((nSpPhase == 1) && (volPhase->phiVarIndex() == 0)) { volPhase->setExistence(VCS_PHASE_EXIST_ALWAYS); } else if (volPhase->totalMoles() > 0.0) { @@ -1021,10 +1011,10 @@ int vcs_Cantera_update_vprob(MultiPhase* mphase, VCS_PROB* vprob) // front of the object. There may be more constraints than there are // elements. But, we know the element abundances are in the front of the // vector. - vprob->set_gai(); + vsolve->set_gai(); // Printout the species information: PhaseID's and mole nums - if (vprob->m_printLvl > 1) { + if (vsolve->m_printLvl > 1) { writeline('=', 80, true, true); writeline('=', 20, false); plogf(" Cantera_to_vprob: START OF PROBLEM STATEMENT "); @@ -1034,15 +1024,15 @@ int vcs_Cantera_update_vprob(MultiPhase* mphase, VCS_PROB* vprob) plogf(" Phase IDs of species\n"); plogf(" species phaseID phaseName "); plogf(" Initial_Estimated_kMols\n"); - for (size_t i = 0; i < vprob->nspecies; i++) { - size_t iphase = vprob->PhaseID[i]; - vcs_VolPhase* VolPhase = vprob->VPhaseList[iphase]; - plogf("%16s %5d %16s", vprob->SpName[i].c_str(), iphase, + for (size_t i = 0; i < vsolve->nspecies; i++) { + size_t iphase = vsolve->PhaseID[i]; + vcs_VolPhase* VolPhase = vsolve->VPhaseList[iphase]; + plogf("%16s %5d %16s", mphase->speciesName(i).c_str(), iphase, VolPhase->PhaseName.c_str()); - if (vprob->SpeciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - plogf(" Volts = %-10.5g\n", vprob->w[i]); + if (vsolve->SpeciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { + plogf(" Volts = %-10.5g\n", vsolve->w[i]); } else { - plogf(" %-10.5g\n", vprob->w[i]); + plogf(" %-10.5g\n", vsolve->w[i]); } } @@ -1052,8 +1042,8 @@ int vcs_Cantera_update_vprob(MultiPhase* mphase, VCS_PROB* vprob) plogf(" PhaseName PhaseNum SingSpec GasPhase EqnState NumSpec"); plogf(" TMolesInert Tmoles(kmol)\n"); - for (size_t iphase = 0; iphase < vprob->NPhase; iphase++) { - vcs_VolPhase* VolPhase = vprob->VPhaseList[iphase]; + for (size_t iphase = 0; iphase < vsolve->NPhase; iphase++) { + vcs_VolPhase* VolPhase = vsolve->VPhaseList[iphase]; plogf("%16s %5d %5d %8d %16s %8d %16e ", VolPhase->PhaseName.c_str(), VolPhase->VP_ID_, VolPhase->m_singleSpecies, VolPhase->m_gasPhase, VolPhase->eos_name(), diff --git a/src/equil/vcs_prep.cpp b/src/equil/vcs_prep.cpp index 1df67050b..9c93bb63c 100644 --- a/src/equil/vcs_prep.cpp +++ b/src/equil/vcs_prep.cpp @@ -7,8 +7,8 @@ // at http://www.cantera.org/license.txt for license and copyright information. #include "cantera/equil/vcs_solve.h" -#include "cantera/equil/vcs_prob.h" #include "cantera/equil/vcs_VolPhase.h" +#include "cantera/equil/MultiPhase.h" namespace Cantera { @@ -195,11 +195,11 @@ int VCS_SOLVE::vcs_prep() return VCS_SUCCESS; } -bool VCS_SOLVE::vcs_wellPosed(VCS_PROB* vprob) +bool VCS_SOLVE::vcs_wellPosed() { double sum = 0.0; - for (size_t e = 0; e < vprob->ne; e++) { - sum += vprob->gai[e]; + for (size_t e = 0; e < ne; e++) { + sum += m_mix->elementMoles(e); } if (sum < 1.0E-20) { plogf("vcs_wellPosed: Element abundance is close to zero\n"); diff --git a/src/equil/vcs_prob.cpp b/src/equil/vcs_prob.cpp index 3a04803f8..3ceb13708 100644 --- a/src/equil/vcs_prob.cpp +++ b/src/equil/vcs_prob.cpp @@ -7,12 +7,13 @@ // This file is part of Cantera. See License.txt in the top-level directory or // at http://www.cantera.org/license.txt for license and copyright information. -#include "cantera/equil/vcs_prob.h" #include "cantera/equil/vcs_VolPhase.h" #include "cantera/equil/vcs_species_thermo.h" #include "cantera/equil/vcs_internal.h" #include "cantera/equil/vcs_defs.h" +#include "cantera/equil/vcs_solve.h" #include "cantera/thermo/MolalityVPSSTP.h" +#include "cantera/equil/MultiPhase.h" #include @@ -21,100 +22,7 @@ using namespace std; namespace Cantera { -VCS_PROB::VCS_PROB(size_t nsp, size_t nel, size_t nph) : - prob_type(VCS_PROBTYPE_TP), - nspecies(nsp), - NSPECIES0(0), - ne(nel), - NE0(0), - NPhase(nph), - NPHASE0(0), - T(298.15), - PresPA(1.0), - Vol(0.0), - // The default is to not expect an initial estimate of the species - // concentrations - iest(-1), - tolmaj(1.0E-8), - tolmin(1.0E-6), - m_Iterations(0), - m_NumBasisOptimizations(0), - m_printLvl(0), - vcs_debug_print_lvl(0) -{ - NSPECIES0 = nspecies; - if (nspecies <= 0) { - throw CanteraError("VCS_PROB::VCS_PROB", - "number of species is zero or neg"); - } - NE0 = ne; - if (ne <= 0) { - throw CanteraError("VCS_PROB::VCS_PROB", - "number of elements is zero or neg"); - } - NPHASE0 = NPhase; - if (NPhase <= 0) { - throw CanteraError("VCS_PROB::VCS_PROB", - "number of phases is zero or neg"); - } - if (nspecies < NPhase) { - throw CanteraError("VCS_PROB::VCS_PROB", - "number of species is less than number of phases"); - } - - m_gibbsSpecies.resize(nspecies, 0.0); - w.resize(nspecies, 0.0); - mf.resize(nspecies, 0.0); - gai.resize(ne, 0.0); - FormulaMatrix.resize(nspecies, ne, 0.0); - SpeciesUnknownType.resize(nspecies, VCS_SPECIES_TYPE_MOLNUM); - VolPM.resize(nspecies, 0.0); - PhaseID.resize(nspecies, npos); - SpName.resize(nspecies, ""); - ElName.resize(ne, ""); - m_elType.resize(ne, VCS_ELEM_TYPE_ABSPOS); - ElActive.resize(ne, 1); - WtSpecies.resize(nspecies, 0.0); - Charge.resize(nspecies, 0.0); - SpeciesThermo.resize(nspecies,0); - for (size_t kspec = 0; kspec < nspecies; kspec++) { - VCS_SPECIES_THERMO* ts_tmp = new VCS_SPECIES_THERMO(0, 0); - if (ts_tmp == 0) { - throw CanteraError("VCS_PROB::VCS_PROB", - "Failed to init a ts struct"); - } - SpeciesThermo[kspec] = ts_tmp; - } - VPhaseList.resize(nph, 0); - for (size_t iphase = 0; iphase < NPhase; iphase++) { - VPhaseList[iphase] = new vcs_VolPhase(); - } -} - -VCS_PROB::~VCS_PROB() -{ - for (size_t i = 0; i < nspecies; i++) { - delete SpeciesThermo[i]; - } - for (size_t iph = 0; iph < NPhase; iph++) { - delete VPhaseList[iph]; - } -} - -void VCS_PROB::resizeElements(size_t nel, int force) -{ - if (force || nel > NE0) { - gai.resize(nel, 0.0); - FormulaMatrix.resize(NSPECIES0, nel, 0.0); - ElName.resize(nel, ""); - m_elType.resize(nel, VCS_ELEM_TYPE_ABSPOS); - ElActive.resize(nel, 1); - NE0 = nel; - ne = std::min(ne, NE0); - } -} - -void VCS_PROB::set_gai() +void VCS_SOLVE::set_gai() { gai.assign(gai.size(), 0.0); for (size_t j = 0; j < ne; j++) { @@ -126,7 +34,7 @@ void VCS_PROB::set_gai() } } -void VCS_PROB::prob_report(int print_lvl) +void VCS_SOLVE::prob_report(int print_lvl) { m_printLvl = print_lvl; @@ -153,7 +61,7 @@ void VCS_PROB::prob_report(int print_lvl) plogf(" Initial_Estimated_Moles Species_Type\n"); for (size_t i = 0; i < nspecies; i++) { vcs_VolPhase* Vphase = VPhaseList[PhaseID[i]]; - plogf("%16s %5d %16s", SpName[i], PhaseID[i], + plogf("%16s %5d %16s", m_mix->speciesName(i), PhaseID[i], Vphase->PhaseName); if (iest >= 0) { plogf(" %-10.5g", w[i]); @@ -194,7 +102,7 @@ void VCS_PROB::prob_report(int print_lvl) plogf(" Target_kmol ElemType ElActive\n"); for (size_t i = 0; i < ne; ++i) { writeline(' ', 26, false); - plogf("%-2.2s", ElName[i]); + plogf("%-2.2s", m_elementName[i]); plogf("%20.12E ", gai[i]); plogf("%3d %3d\n", m_elType[i], ElActive[i]); } @@ -207,7 +115,7 @@ void VCS_PROB::prob_report(int print_lvl) Vphase->setState_TP(T, PresPA); for (size_t kindex = 0; kindex < Vphase->nSpecies(); kindex++) { size_t kglob = Vphase->spGlobalIndexVCS(kindex); - plogf("%16s ", SpName[kglob]); + plogf("%16s ", m_mix->speciesName(kglob)); if (kindex == 0) { plogf("%16s", Vphase->PhaseName); } else { @@ -227,10 +135,9 @@ void VCS_PROB::prob_report(int print_lvl) } } -void VCS_PROB::addPhaseElements(vcs_VolPhase* volPhase) +void VCS_SOLVE::addPhaseElements(vcs_VolPhase* volPhase) { size_t neVP = volPhase->nElemConstraints(); - // Loop through the elements in the vol phase object for (size_t eVP = 0; eVP < neVP; eVP++) { size_t foundPos = npos; @@ -239,7 +146,7 @@ void VCS_PROB::addPhaseElements(vcs_VolPhase* volPhase) // Search for matches with the existing elements. If found, then fill in // the entry in the global mapping array. for (size_t e = 0; e < ne; e++) { - std::string en = ElName[e]; + std::string en = m_elementName[e]; if (!strcmp(enVP.c_str(), en.c_str())) { volPhase->setElemGlobalIndex(eVP, e); foundPos = e; @@ -254,22 +161,32 @@ void VCS_PROB::addPhaseElements(vcs_VolPhase* volPhase) } } -size_t VCS_PROB::addElement(const char* elNameNew, int elType, int elactive) +size_t VCS_SOLVE::addElement(const char* elNameNew, int elType, int elactive) { if (!elNameNew) { - throw CanteraError("VCS_PROB::addElement", + throw CanteraError("VCS_SOLVE::addElement", "error: element must have a name"); } - size_t nel = ne + 1; - resizeElements(nel, 1); - ne = nel; - ElName[ne-1] = elNameNew; - m_elType[ne-1] = elType; - ElActive[ne-1] = elactive; + ne++; + m_numElemConstraints++; + m_numComponents++; + + gai.push_back(0.0); + FormulaMatrix.resize(NSPECIES0, ne, 0.0); + m_formulaMatrix.resize(NSPECIES0, ne); + m_stoichCoeffRxnMatrix.resize(ne, NSPECIES0, 0.0); + m_elType.push_back(elType); + ElActive.push_back(elactive); + m_elementActive.push_back(elactive); + m_elemAbundances.push_back(0.0); + m_elemAbundancesGoal.push_back(0.0); + m_elementMapIndex.push_back(0); + m_elementName.push_back(elNameNew); + NE0 = ne; return ne - 1; } -size_t VCS_PROB::addOnePhaseSpecies(vcs_VolPhase* volPhase, size_t k, size_t kT) +size_t VCS_SOLVE::addOnePhaseSpecies(vcs_VolPhase* volPhase, size_t k, size_t kT) { if (kT > nspecies) { // Need to expand the number of species here @@ -289,9 +206,4 @@ size_t VCS_PROB::addOnePhaseSpecies(vcs_VolPhase* volPhase, size_t k, size_t kT) return kT; } -void VCS_PROB::setDebugPrintLvl(int lvl) -{ - vcs_debug_print_lvl = lvl; -} - } diff --git a/src/equil/vcs_solve.cpp b/src/equil/vcs_solve.cpp index 99a0758a7..023425194 100644 --- a/src/equil/vcs_solve.cpp +++ b/src/equil/vcs_solve.cpp @@ -9,10 +9,10 @@ #include "cantera/equil/vcs_solve.h" #include "cantera/base/ctexceptions.h" #include "cantera/base/stringUtils.h" -#include "cantera/equil/vcs_prob.h" #include "cantera/equil/vcs_VolPhase.h" #include "cantera/equil/vcs_species_thermo.h" #include "cantera/base/clockWC.h" +#include "cantera/equil/MultiPhase.h" using namespace std; @@ -21,10 +21,27 @@ namespace Cantera int vcs_timing_print_lvl = 1; -VCS_SOLVE::VCS_SOLVE() : - NSPECIES0(0), - NPHASE0(0), - m_numSpeciesTot(0), +VCS_SOLVE::VCS_SOLVE(size_t nspecies_, size_t nelements_, size_t nphase_) : + prob_type(VCS_PROBTYPE_TP), + nspecies(nspecies_), + ne(0), + NE0(0), + NPhase(nphase_), + T(298.15), + PresPA(1.0), + Vol(0.0), + // The default is to not expect an initial estimate of the species + // concentrations + iest(-1), + tolmaj(1.0E-8), + tolmin(1.0E-6), + m_Iterations(0), + m_NumBasisOptimizations(0), + m_printLvl(0), + vcs_debug_print_lvl(0), + NSPECIES0(nspecies_), + NPHASE0(nphase_), + m_numSpeciesTot(nspecies_), m_numElemConstraints(0), m_numComponents(0), m_numRxnTot(0), @@ -49,37 +66,30 @@ VCS_SOLVE::VCS_SOLVE() : m_debug_print_lvl(0), m_timing_print_lvl(1) { -} - -void VCS_SOLVE::vcs_initSizes(const size_t nspecies0, const size_t nelements, - const size_t nphase0) -{ - if (NSPECIES0 != 0) { - if ((nspecies0 != NSPECIES0) || (nelements != m_numElemConstraints) || (nphase0 != NPHASE0)) { - vcs_delete_memory(); - } else { - return; - } + m_gibbsSpecies.resize(nspecies, 0.0); + w.resize(nspecies, 0.0); + mf.resize(nspecies, 0.0); + SpeciesUnknownType.resize(nspecies, VCS_SPECIES_TYPE_MOLNUM); + VolPM.resize(nspecies, 0.0); + PhaseID.resize(nspecies, npos); + WtSpecies.resize(nspecies, 0.0); + Charge.resize(nspecies, 0.0); + SpeciesThermo.resize(nspecies,0); + for (size_t kspec = 0; kspec < nspecies; kspec++) { + SpeciesThermo[kspec] = new VCS_SPECIES_THERMO(0, 0); + } + VPhaseList.resize(NPhase, 0); + for (size_t iphase = 0; iphase < NPhase; iphase++) { + VPhaseList[iphase] = new vcs_VolPhase(); } - NSPECIES0 = nspecies0; - NPHASE0 = nphase0; - m_numSpeciesTot = nspecies0; - m_numElemConstraints = nelements; - m_numComponents = nelements; - string ser = "VCS_SOLVE: ERROR:\n\t"; - if (nspecies0 <= 0) { + if (nspecies_ <= 0) { plogf("%s Number of species is nonpositive\n", ser); throw CanteraError("VCS_SOLVE()", ser + " Number of species is nonpositive\n"); } - if (nelements <= 0) { - plogf("%s Number of elements is nonpositive\n", ser); - throw CanteraError("VCS_SOLVE()", ser + - " Number of species is nonpositive\n"); - } - if (nphase0 <= 0) { + if (nphase_ <= 0) { plogf("%s Number of phases is nonpositive\n", ser); throw CanteraError("VCS_SOLVE()", ser + " Number of species is nonpositive\n"); @@ -89,79 +99,68 @@ void VCS_SOLVE::vcs_initSizes(const size_t nspecies0, const size_t nelements, * We will initialize sc[] to note the fact that it needs to be * filled with meaningful information. */ - m_stoichCoeffRxnMatrix.resize(nelements, nspecies0, 0.0); - m_scSize.resize(nspecies0, 0.0); - m_spSize.resize(nspecies0, 1.0); - m_SSfeSpecies.resize(nspecies0, 0.0); - m_feSpecies_new.resize(nspecies0, 0.0); - m_molNumSpecies_old.resize(nspecies0, 0.0); - m_speciesUnknownType.resize(nspecies0, VCS_SPECIES_TYPE_MOLNUM); - m_deltaMolNumPhase.resize(nphase0, nspecies0, 0.0); - m_phaseParticipation.resize(nphase0, nspecies0, 0); - m_phasePhi.resize(nphase0, 0.0); - m_molNumSpecies_new.resize(nspecies0, 0.0); - m_deltaGRxn_new.resize(nspecies0, 0.0); - m_deltaGRxn_old.resize(nspecies0, 0.0); - m_deltaGRxn_Deficient.resize(nspecies0, 0.0); - m_deltaGRxn_tmp.resize(nspecies0, 0.0); - m_deltaMolNumSpecies.resize(nspecies0, 0.0); - m_feSpecies_old.resize(nspecies0, 0.0); - m_elemAbundances.resize(nelements, 0.0); - m_elemAbundancesGoal.resize(nelements, 0.0); - m_tPhaseMoles_old.resize(nphase0, 0.0); - m_tPhaseMoles_new.resize(nphase0, 0.0); - m_deltaPhaseMoles.resize(nphase0, 0.0); - m_TmpPhase.resize(nphase0, 0.0); - m_TmpPhase2.resize(nphase0, 0.0); - m_formulaMatrix.resize(nspecies0, nelements); - TPhInertMoles.resize(nphase0, 0.0); + m_scSize.resize(NSPECIES0, 0.0); + m_spSize.resize(NSPECIES0, 1.0); + m_SSfeSpecies.resize(NSPECIES0, 0.0); + m_feSpecies_new.resize(NSPECIES0, 0.0); + m_molNumSpecies_old.resize(NSPECIES0, 0.0); + m_speciesUnknownType.resize(NSPECIES0, VCS_SPECIES_TYPE_MOLNUM); + m_deltaMolNumPhase.resize(NPHASE0, NSPECIES0, 0.0); + m_phaseParticipation.resize(NPHASE0, NSPECIES0, 0); + m_phasePhi.resize(NPHASE0, 0.0); + m_molNumSpecies_new.resize(NSPECIES0, 0.0); + m_deltaGRxn_new.resize(NSPECIES0, 0.0); + m_deltaGRxn_old.resize(NSPECIES0, 0.0); + m_deltaGRxn_Deficient.resize(NSPECIES0, 0.0); + m_deltaGRxn_tmp.resize(NSPECIES0, 0.0); + m_deltaMolNumSpecies.resize(NSPECIES0, 0.0); + m_feSpecies_old.resize(NSPECIES0, 0.0); + m_tPhaseMoles_old.resize(NPHASE0, 0.0); + m_tPhaseMoles_new.resize(NPHASE0, 0.0); + m_deltaPhaseMoles.resize(NPHASE0, 0.0); + m_TmpPhase.resize(NPHASE0, 0.0); + m_TmpPhase2.resize(NPHASE0, 0.0); + TPhInertMoles.resize(NPHASE0, 0.0); // ind[] is an index variable that keep track of solution vector rotations. - m_speciesMapIndex.resize(nspecies0, 0); - m_speciesLocalPhaseIndex.resize(nspecies0, 0); - - // IndEl[] is an index variable that keep track of element vector rotations. - m_elementMapIndex.resize(nelements, 0); + m_speciesMapIndex.resize(NSPECIES0, 0); + m_speciesLocalPhaseIndex.resize(NSPECIES0, 0); // ir[] is an index vector that keeps track of the irxn to species mapping. // We can't fill it in until we know the number of c components in the // problem - m_indexRxnToSpecies.resize(nspecies0, 0); + m_indexRxnToSpecies.resize(NSPECIES0, 0); // Initialize all species to be major species - m_speciesStatus.resize(nspecies0, 1); + m_speciesStatus.resize(NSPECIES0, 1); - m_SSPhase.resize(2*nspecies0, 0); - m_phaseID.resize(nspecies0, 0); - m_numElemConstraints = nelements; - m_elementName.resize(nelements); - m_speciesName.resize(nspecies0); - m_elType.resize(nelements, VCS_ELEM_TYPE_ABSPOS); - m_elementActive.resize(nelements, 1); + m_SSPhase.resize(2*NSPECIES0, 0); + m_phaseID.resize(NSPECIES0, 0); + m_speciesName.resize(NSPECIES0); // space for activity coefficients for all species. Set it equal to one. - m_actConventionSpecies.resize(nspecies0, 0); - m_phaseActConvention.resize(nphase0, 0); - m_lnMnaughtSpecies.resize(nspecies0, 0.0); - m_actCoeffSpecies_new.resize(nspecies0, 1.0); - m_actCoeffSpecies_old.resize(nspecies0, 1.0); - m_wtSpecies.resize(nspecies0, 0.0); - m_chargeSpecies.resize(nspecies0, 0.0); - m_speciesThermoList.resize(nspecies0, (VCS_SPECIES_THERMO*)0); + m_actConventionSpecies.resize(NSPECIES0, 0); + m_phaseActConvention.resize(NPHASE0, 0); + m_lnMnaughtSpecies.resize(NSPECIES0, 0.0); + m_actCoeffSpecies_new.resize(NSPECIES0, 1.0); + m_actCoeffSpecies_old.resize(NSPECIES0, 1.0); + m_wtSpecies.resize(NSPECIES0, 0.0); + m_chargeSpecies.resize(NSPECIES0, 0.0); + m_speciesThermoList.resize(NSPECIES0, (VCS_SPECIES_THERMO*)0); // Phase Info - m_VolPhaseList.resize(nphase0, 0); - for (size_t iph = 0; iph < nphase0; iph++) { + m_VolPhaseList.resize(NPHASE0, 0); + for (size_t iph = 0; iph < NPHASE0; iph++) { m_VolPhaseList[iph] = new vcs_VolPhase(this); } // For Future expansion m_useActCoeffJac = true; if (m_useActCoeffJac) { - m_np_dLnActCoeffdMolNum.resize(nspecies0, nspecies0, 0.0); + m_np_dLnActCoeffdMolNum.resize(NSPECIES0, NSPECIES0, 0.0); } - m_PMVolumeSpecies.resize(nspecies0, 0.0); + m_PMVolumeSpecies.resize(NSPECIES0, 0.0); // counters kept within vcs m_VCount = new VCS_COUNTERS(); @@ -181,18 +180,25 @@ VCS_SOLVE::~VCS_SOLVE() void VCS_SOLVE::vcs_delete_memory() { - size_t nspecies = m_numSpeciesTot; - - for (size_t j = 0; j < m_numPhases; j++) { + for (size_t j = 0; j < m_VolPhaseList.size(); j++) { delete m_VolPhaseList[j]; m_VolPhaseList[j] = 0; } - for (size_t j = 0; j < nspecies; j++) { + for (size_t j = 0; j < m_speciesThermoList.size(); j++) { delete m_speciesThermoList[j]; m_speciesThermoList[j] = 0; } + for (size_t i = 0; i < SpeciesThermo.size(); i++) { + delete SpeciesThermo[i]; + SpeciesThermo[i] = 0; + } + for (size_t iph = 0; iph < VPhaseList.size(); iph++) { + delete VPhaseList[iph]; + VPhaseList[iph] = 0; + } + delete m_VCount; m_VCount = 0; @@ -200,15 +206,18 @@ void VCS_SOLVE::vcs_delete_memory() NPHASE0 = 0; m_numElemConstraints = 0; m_numComponents = 0; + } -int VCS_SOLVE::vcs(VCS_PROB* vprob, int ifunc, int ipr, int ip1, int maxit) +int VCS_SOLVE::vcs(int ifunc, int ipr, int ip1, int maxit) { int retn = 0, iconv = 0; clockWC tickTock; int iprintTime = std::max(ipr, ip1); - iprintTime = std::min(iprintTime, m_timing_print_lvl); + //iprintTime = std::min(iprintTime, m_timing_print_lvl); + // writelog("iprintTime = {}; m_timing_print_lvl = {}; ipr = {}; ip1 = {}\n", + // iprintTime, m_timing_print_lvl, ipr, ip1); if (ifunc > 2) { plogf("vcs: Unrecognized value of ifunc, %d: bailing!\n", @@ -217,21 +226,9 @@ int VCS_SOLVE::vcs(VCS_PROB* vprob, int ifunc, int ipr, int ip1, int maxit) } if (ifunc == 0) { - // This function is called to create the private data using the public - // data. - size_t nspecies0 = vprob->nspecies + 10; - size_t nelements0 = vprob->ne; - size_t nphase0 = vprob->NPhase; - vcs_initSizes(nspecies0, nelements0, nphase0); - - if (retn != 0) { - plogf("vcs_priv_alloc returned a bad status, %d: bailing!\n", - retn); - return retn; - } // This function is called to copy the public data and the current // problem specification into the current object's data structure. - retn = vcs_prob_specifyFully(vprob); + retn = vcs_prob_specifyFully(); if (retn != 0) { plogf("vcs_pub_to_priv returned a bad status, %d: bailing!\n", retn); @@ -250,7 +247,7 @@ int VCS_SOLVE::vcs(VCS_PROB* vprob, int ifunc, int ipr, int ip1, int maxit) if (ifunc == 1) { // This function is called to copy the current problem into the current // object's data structure. - retn = vcs_prob_specify(vprob); + retn = vcs_prob_specify(); if (retn != 0) { plogf("vcs_prob_specify returned a bad status, %d: bailing!\n", retn); @@ -267,7 +264,7 @@ int VCS_SOLVE::vcs(VCS_PROB* vprob, int ifunc, int ipr, int ip1, int maxit) } // Check to see if the current problem is well posed. - if (!vcs_wellPosed(vprob)) { + if (!vcs_wellPosed()) { plogf("vcs has determined the problem is not well posed: Bailing\n"); return VCS_PUB_BAD; } @@ -279,16 +276,14 @@ int VCS_SOLVE::vcs(VCS_PROB* vprob, int ifunc, int ipr, int ip1, int maxit) // problem types will go in at this level. For example, solving for // fixed T, V problems will involve a 2x2 Newton's method, using loops // over vcs_TP() to calculate the residual and Jacobian) - iconv = vcs_TP(ipr, ip1, maxit, vprob->T, vprob->PresPA); + iconv = vcs_TP(ipr, ip1, maxit, T, PresPA); // If requested to print anything out, go ahead and do so; if (ipr > 0) { vcs_report(iconv); } - // Copy the results of the run back to the VCS_PROB structure, which is - // returned to the user. - vcs_prob_update(vprob); + vcs_prob_update(); } // Report on the time if requested to do so @@ -308,45 +303,34 @@ int VCS_SOLVE::vcs(VCS_PROB* vprob, int ifunc, int ipr, int ip1, int maxit) return iconv; } -int VCS_SOLVE::vcs_prob_specifyFully(const VCS_PROB* pub) +int VCS_SOLVE::vcs_prob_specifyFully() { const char* ser = "vcs_pub_to_priv ERROR :ill defined interface -> bailout:\n\t"; // First Check to see whether we have room for the current problem size - size_t nspecies = pub->nspecies; if (NSPECIES0 < nspecies) { plogf("%sPrivate Data is dimensioned too small\n", ser); return VCS_PUB_BAD; } - size_t nph = pub->NPhase; + size_t nph = NPhase; if (NPHASE0 < nph) { plogf("%sPrivate Data is dimensioned too small\n", ser); return VCS_PUB_BAD; } - size_t nelements = pub->ne; - if (m_numElemConstraints < nelements) { - plogf("%sPrivate Data is dimensioned too small\n", ser); - return VCS_PUB_BAD; - } // OK, We have room. Now, transfer the integer numbers - m_numElemConstraints = nelements; m_numSpeciesTot = nspecies; m_numSpeciesRdc = m_numSpeciesTot; - // nc = number of components -> will be determined later. but set it to its - // maximum possible value here. - m_numComponents = nelements; - // m_numRxnTot = number of noncomponents, also equal to the number of // reactions. Note, it's possible that the number of elements is greater // than the number of species. In that case set the number of reactions to // zero. - if (nelements > nspecies) { + if (ne > nspecies) { m_numRxnTot = 0; } else { - m_numRxnTot = nspecies - nelements; + m_numRxnTot = nspecies - ne; } m_numRxnRdc = m_numRxnTot; @@ -356,34 +340,34 @@ int VCS_SOLVE::vcs_prob_specifyFully(const VCS_PROB* pub) // NPhase = number of phases m_numPhases = nph; - m_debug_print_lvl = pub->vcs_debug_print_lvl; + m_debug_print_lvl = vcs_debug_print_lvl; // FormulaMatrix[] -> Copy the formula matrix over for (size_t i = 0; i < nspecies; i++) { bool nonzero = false; - for (size_t j = 0; j < nelements; j++) { - if (pub->FormulaMatrix(i,j) != 0.0) { + for (size_t j = 0; j < ne; j++) { + if (FormulaMatrix(i,j) != 0.0) { nonzero = true; } - m_formulaMatrix(i,j) = pub->FormulaMatrix(i,j); + m_formulaMatrix(i,j) = FormulaMatrix(i,j); } if (!nonzero) { plogf("vcs_prob_specifyFully:: species %d %s has a zero formula matrix!\n", i, - pub->SpName[i]); + m_speciesName[i]); return VCS_PUB_BAD; } } // Copy over the species molecular weights - m_wtSpecies = pub->WtSpecies; + m_wtSpecies = WtSpecies; // Copy over the charges - m_chargeSpecies = pub->Charge; + m_chargeSpecies = Charge; // Copy the VCS_SPECIES_THERMO structures for (size_t kspec = 0; kspec < nspecies; kspec++) { delete m_speciesThermoList[kspec]; - VCS_SPECIES_THERMO* spf = pub->SpeciesThermo[kspec]; + VCS_SPECIES_THERMO* spf = SpeciesThermo[kspec]; m_speciesThermoList[kspec] = spf->duplMyselfAsVCS_SPECIES_THERMO(); if (m_speciesThermoList[kspec] == NULL) { plogf(" duplMyselfAsVCS_SPECIES_THERMO returned an error!\n"); @@ -392,31 +376,31 @@ int VCS_SOLVE::vcs_prob_specifyFully(const VCS_PROB* pub) } // Copy the species unknown type - m_speciesUnknownType = pub->SpeciesUnknownType; + m_speciesUnknownType = SpeciesUnknownType; // iest => Do we have an initial estimate of the species mole numbers ? - m_doEstimateEquil = pub->iest; + m_doEstimateEquil = iest; // w[] -> Copy the equilibrium mole number estimate if it exists. - if (pub->w.size() != 0) { - m_molNumSpecies_old = pub->w; + if (w.size() != 0) { + m_molNumSpecies_old = w; } else { m_doEstimateEquil = -1; m_molNumSpecies_old.assign(m_molNumSpecies_old.size(), 0.0); } // Formulate the Goal Element Abundance Vector - if (pub->gai.size() != 0) { - for (size_t i = 0; i < nelements; i++) { - m_elemAbundancesGoal[i] = pub->gai[i]; - if (pub->m_elType[i] == VCS_ELEM_TYPE_LATTICERATIO && m_elemAbundancesGoal[i] < 1.0E-10) { + if (gai.size() != 0) { + for (size_t i = 0; i < ne; i++) { + m_elemAbundancesGoal[i] = gai[i]; + if (m_elType[i] == VCS_ELEM_TYPE_LATTICERATIO && m_elemAbundancesGoal[i] < 1.0E-10) { m_elemAbundancesGoal[i] = 0.0; } } } else { if (m_doEstimateEquil == 0) { double sum = 0; - for (size_t j = 0; j < nelements; j++) { + for (size_t j = 0; j < ne; j++) { m_elemAbundancesGoal[j] = 0.0; for (size_t kspec = 0; kspec < nspecies; kspec++) { if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { @@ -424,7 +408,7 @@ int VCS_SOLVE::vcs_prob_specifyFully(const VCS_PROB* pub) m_elemAbundancesGoal[j] += m_formulaMatrix(kspec,j) * m_molNumSpecies_old[kspec]; } } - if (pub->m_elType[j] == VCS_ELEM_TYPE_LATTICERATIO && m_elemAbundancesGoal[j] < 1.0E-10 * sum) { + if (m_elType[j] == VCS_ELEM_TYPE_LATTICERATIO && m_elemAbundancesGoal[j] < 1.0E-10 * sum) { m_elemAbundancesGoal[j] = 0.0; } } @@ -441,26 +425,26 @@ int VCS_SOLVE::vcs_prob_specifyFully(const VCS_PROB* pub) // DelTPhMoles[] // // T, Pres, copy over here - if (pub->T > 0.0) { - m_temperature = pub->T; + if (T > 0.0) { + m_temperature = T; } else { m_temperature = 293.15; } - if (pub->PresPA > 0.0) { - m_pressurePA = pub->PresPA; + if (PresPA > 0.0) { + m_pressurePA = PresPA; } else { m_pressurePA = OneAtm; } // TPhInertMoles[] -> must be copied over here for (size_t iph = 0; iph < nph; iph++) { - vcs_VolPhase* Vphase = pub->VPhaseList[iph]; + vcs_VolPhase* Vphase = VPhaseList[iph]; TPhInertMoles[iph] = Vphase->totalMolesInert(); } // tolerance requirements -> copy them over here and later - m_tolmaj = pub->tolmaj; - m_tolmin = pub->tolmin; + m_tolmaj = tolmaj; + m_tolmin = tolmin; m_tolmaj2 = 0.01 * m_tolmaj; m_tolmin2 = 0.01 * m_tolmin; @@ -471,7 +455,7 @@ int VCS_SOLVE::vcs_prob_specifyFully(const VCS_PROB* pub) } // IndEl[] is an index variable that keep track of element vector rotations. - for (size_t i = 0; i < nelements; i++) { + for (size_t i = 0; i < ne; i++) { m_elementMapIndex[i] = i; } @@ -482,10 +466,10 @@ int VCS_SOLVE::vcs_prob_specifyFully(const VCS_PROB* pub) // PhaseID: Fill in the species to phase mapping. Check for bad values at // the same time. - if (pub->PhaseID.size() != 0) { + if (PhaseID.size() != 0) { std::vector numPhSp(nph, 0); for (size_t kspec = 0; kspec < nspecies; kspec++) { - size_t iph = pub->PhaseID[kspec]; + size_t iph = PhaseID[kspec]; if (iph >= nph) { plogf("%sSpecies to Phase Mapping, PhaseID, has a bad value\n", ser); @@ -493,15 +477,15 @@ int VCS_SOLVE::vcs_prob_specifyFully(const VCS_PROB* pub) plogf("\tAllowed values: 0 to %d\n", nph - 1); return VCS_PUB_BAD; } - m_phaseID[kspec] = pub->PhaseID[kspec]; + m_phaseID[kspec] = PhaseID[kspec]; m_speciesLocalPhaseIndex[kspec] = numPhSp[iph]; numPhSp[iph]++; } for (size_t iph = 0; iph < nph; iph++) { - vcs_VolPhase* Vphase = pub->VPhaseList[iph]; + vcs_VolPhase* Vphase = VPhaseList[iph]; if (numPhSp[iph] != Vphase->nSpecies()) { - plogf("%sNumber of species in phase %d, %s, doesn't match\n", - ser, iph, Vphase->PhaseName); + plogf("%sNumber of species in phase %d, %s, doesn't match (%d != %d) [vphase = %d]\n", + ser, iph, Vphase->PhaseName, numPhSp[iph], Vphase->nSpecies(), (size_t) Vphase); return VCS_PUB_BAD; } } @@ -517,25 +501,7 @@ int VCS_SOLVE::vcs_prob_specifyFully(const VCS_PROB* pub) } } - // Copy over the element types - m_elType.resize(nelements, VCS_ELEM_TYPE_ABSPOS); - m_elementActive.resize(nelements, 1); - - // Copy over the element names and types - for (size_t i = 0; i < nelements; i++) { - m_elementName[i] = pub->ElName[i]; - m_elType[i] = pub->m_elType[i]; - m_elementActive[i] = pub->ElActive[i]; - if (!strncmp(m_elementName[i].c_str(), "cn_", 3)) { - m_elType[i] = VCS_ELEM_TYPE_CHARGENEUTRALITY; - if (pub->m_elType[i] != VCS_ELEM_TYPE_CHARGENEUTRALITY) { - throw CanteraError("VCS_SOLVE::vcs_prob_specifyFully", - "we have an inconsistency!"); - } - } - } - - for (size_t i = 0; i < nelements; i++) { + for (size_t i = 0; i < ne; i++) { if (m_elType[i] == VCS_ELEM_TYPE_CHARGENEUTRALITY) { if (m_elemAbundancesGoal[i] != 0.0) { if (fabs(m_elemAbundancesGoal[i]) > 1.0E-9) { @@ -557,13 +523,13 @@ int VCS_SOLVE::vcs_prob_specifyFully(const VCS_PROB* pub) // Copy over the species names for (size_t i = 0; i < nspecies; i++) { - m_speciesName[i] = pub->SpName[i]; + m_speciesName[i] = m_mix->speciesName(i); } // Copy over all of the phase information. Use the object's assignment // operator for (size_t iph = 0; iph < nph; iph++) { - *m_VolPhaseList[iph] = *pub->VPhaseList[iph]; + *m_VolPhaseList[iph] = *VPhaseList[iph]; // Fix up the species thermo pointer in the vcs_SpeciesThermo object. It // should point to the species thermo pointer in the private data space. @@ -597,56 +563,56 @@ int VCS_SOLVE::vcs_prob_specifyFully(const VCS_PROB* pub) } // Copy the title info - if (pub->Title.size() == 0) { + if (Title.size() == 0) { m_title = "Unspecified Problem Title"; } else { - m_title = pub->Title; + m_title = Title; } // Copy the volume info - m_totalVol = pub->Vol; + m_totalVol = Vol; if (m_PMVolumeSpecies.size() != 0) { - m_PMVolumeSpecies = pub->VolPM; + m_PMVolumeSpecies = VolPM; } // Return the success flag return VCS_SUCCESS; } -int VCS_SOLVE::vcs_prob_specify(const VCS_PROB* pub) +int VCS_SOLVE::vcs_prob_specify() { string yo("vcs_prob_specify ERROR: "); int retn = VCS_SUCCESS; - m_temperature = pub->T; - m_pressurePA = pub->PresPA; - m_doEstimateEquil = pub->iest; - m_totalVol = pub->Vol; - m_tolmaj = pub->tolmaj; - m_tolmin = pub->tolmin; + m_temperature = T; + m_pressurePA = PresPA; + m_doEstimateEquil = iest; + m_totalVol = Vol; + m_tolmaj = tolmaj; + m_tolmin = tolmin; m_tolmaj2 = 0.01 * m_tolmaj; m_tolmin2 = 0.01 * m_tolmin; for (size_t kspec = 0; kspec < m_numSpeciesTot; ++kspec) { size_t k = m_speciesMapIndex[kspec]; - m_molNumSpecies_old[kspec] = pub->w[k]; - m_molNumSpecies_new[kspec] = pub->mf[k]; - m_feSpecies_old[kspec] = pub->m_gibbsSpecies[k]; + m_molNumSpecies_old[kspec] = w[k]; + m_molNumSpecies_new[kspec] = mf[k]; + m_feSpecies_old[kspec] = m_gibbsSpecies[k]; } // Transfer the element abundance goals to the solve object for (size_t i = 0; i < m_numElemConstraints; i++) { size_t j = m_elementMapIndex[i]; - m_elemAbundancesGoal[i] = pub->gai[j]; + m_elemAbundancesGoal[i] = gai[j]; } // Try to do the best job at guessing at the title - if (pub->Title.size() == 0) { + if (Title.size() == 0) { if (m_title.size() == 0) { m_title = "Unspecified Problem Title"; } } else { - m_title = pub->Title; + m_title = Title; } // Copy over the phase information. For each entry in the phase structure, @@ -655,7 +621,7 @@ int VCS_SOLVE::vcs_prob_specify(const VCS_PROB* pub) bool status_change = false; for (size_t iph = 0; iph < m_numPhases; iph++) { vcs_VolPhase* vPhase = m_VolPhaseList[iph]; - vcs_VolPhase* pub_phase_ptr = pub->VPhaseList[iph]; + vcs_VolPhase* pub_phase_ptr = VPhaseList[iph]; if (vPhase->VP_ID_ != pub_phase_ptr->VP_ID_) { plogf("%sPhase numbers have changed:%d %d\n", @@ -715,7 +681,7 @@ int VCS_SOLVE::vcs_prob_specify(const VCS_PROB* pub) return retn; } -int VCS_SOLVE::vcs_prob_update(VCS_PROB* pub) +int VCS_SOLVE::vcs_prob_update() { size_t k1 = 0; vcs_tmoles(); @@ -733,21 +699,21 @@ int VCS_SOLVE::vcs_prob_update(VCS_PROB* pub) } // Switch the species data back from K1 into I - if (pub->SpeciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - pub->w[i] = m_molNumSpecies_old[k1]; + if (SpeciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { + w[i] = m_molNumSpecies_old[k1]; } else { - pub->w[i] = 0.0; + w[i] = 0.0; } - pub->m_gibbsSpecies[i] = m_feSpecies_old[k1]; - pub->VolPM[i] = m_PMVolumeSpecies[k1]; + m_gibbsSpecies[i] = m_feSpecies_old[k1]; + VolPM[i] = m_PMVolumeSpecies[k1]; } - pub->T = m_temperature; - pub->PresPA = m_pressurePA; - pub->Vol = m_totalVol; + T = m_temperature; + PresPA = m_pressurePA; + Vol = m_totalVol; size_t kT = 0; - for (size_t iph = 0; iph < pub->NPhase; iph++) { - vcs_VolPhase* pubPhase = pub->VPhaseList[iph]; + for (size_t iph = 0; iph < NPhase; iph++) { + vcs_VolPhase* pubPhase = VPhaseList[iph]; vcs_VolPhase* vPhase = m_VolPhaseList[iph]; pubPhase->setTotalMolesInert(vPhase->totalMolesInert()); pubPhase->setTotalMoles(vPhase->totalMoles()); @@ -759,7 +725,7 @@ int VCS_SOLVE::vcs_prob_update(VCS_PROB* pub) const vector_fp & mfVector = pubPhase->moleFractions(); for (size_t k = 0; k < pubPhase->nSpecies(); k++) { kT = pubPhase->spGlobalIndexVCS(k); - pub->mf[kT] = mfVector[k]; + mf[kT] = mfVector[k]; if (pubPhase->phiVarIndex() == k) { k1 = vPhase->spGlobalIndexVCS(k); double tmp = m_molNumSpecies_old[k1]; @@ -770,13 +736,13 @@ int VCS_SOLVE::vcs_prob_update(VCS_PROB* pub) } } - if (! vcs_doubleEqual(pub->mf[kT], vPhase->molefraction(k))) { + if (! vcs_doubleEqual(mf[kT], vPhase->molefraction(k))) { throw CanteraError("VCS_SOLVE::vcs_prob_update", "We have an inconsistency in mole fraction, {} {}", - pub->mf[kT], vPhase->molefraction(k)); + mf[kT], vPhase->molefraction(k)); } if (pubPhase->speciesUnknownType(k) != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - sumMoles += pub->w[kT]; + sumMoles += w[kT]; } } if (! vcs_doubleEqual(sumMoles, vPhase->totalMoles())) { @@ -786,8 +752,8 @@ int VCS_SOLVE::vcs_prob_update(VCS_PROB* pub) } } - pub->m_Iterations = m_VCount->Its; - pub->m_NumBasisOptimizations = m_VCount->Basis_Opts; + m_Iterations = m_VCount->Its; + m_NumBasisOptimizations = m_VCount->Basis_Opts; return VCS_SUCCESS; }