research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

[Equil] Refactor to eliminate class VCS_PROB

Browse source 
Move data in the VCS_PROB class to VCS_SOLVE
This commit is contained in:
Ray Speth 2017-08-16 18:52:36 -04:00
parent da801f4cbc
commit 8522095dea
9 changed files with 550 additions and 787 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

41
include/cantera/equil/vcs_MultiPhaseEquil.h
View file

@ -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;
};

3
include/cantera/equil/vcs_VolPhase.h
View file

@ -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<size_t> m_elemGlobalIndex;
//! Number of species in the phase

1
include/cantera/equil/vcs_internal.h
View file

@ -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
/*!

269
include/cantera/equil/vcs_prob.h
View file

@ -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<size_t> PhaseID;
//! Vector of strings containing the species names
std::vector<std::string> SpName;
//! vector of strings containing the element names
std::vector<std::string> 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<vcs_VolPhase*> 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<VCS_SPECIES_THERMO*> 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

253
include/cantera/equil/vcs_solve.h
View file

@ -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<size_t> 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<vcs_VolPhase*> 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<VCS_SPECIES_THERMO*> 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<size_t> 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<size_t> m_elementMapIndex;

234
src/equil/vcs_MultiPhaseEquil.cpp
View file

@ -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(),

8
src/equil/vcs_prep.cpp
View file

@ -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");

144
src/equil/vcs_prob.cpp
View file

@ -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 <cstdio>
@ -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;
}
}

384
src/equil/vcs_solve.cpp
View file

@ -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<size_t> 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;
}