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

[Equil] Standardize VCS solver to always use MKS units

Browse source
This commit is contained in:
Ray Speth 2016-04-03 21:04:09 -04:00
parent a60217cfc6
commit 511a39fc73
16 changed files with 63 additions and 417 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

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

@ -292,9 +292,6 @@ public:
const ThermoPhase* ptrThermoPhase() const;
//! Return the total moles in the phase
/*!
* Units -> depends on VCS_UnitsFormat variable. Cantera -> J/kmol
*/
double totalMoles() const;
//! Returns the mole fraction of the kspec species
@ -586,31 +583,6 @@ public:
*/
size_t ChargeNeutralityElement;
//! Units for the chemical potential data, pressure data, volume,
//! and species amounts
/*!
* All internally stored quantities will have these units. Also, printed
* quantities will display in these units. Input quantities are expected in
* these units.
*
* | | | Chem_Pot | Pres | vol | moles|
* |---|--------------------|-------------------------|------|------|------|
* |-1 | VCS_UNITS_KCALMOL | kcal/gmol | Pa | m**3 | kmol |
* | 0 | VCS_UNITS_UNITLESS | MU / RT -> no units | Pa | m**3 | kmol |
* | 1 | VCS_UNITS_KJMOL | kJ / gmol | Pa | m**3 | kmol |
* | 2 | VCS_UNITS_KELVIN | KELVIN -> MU / R | Pa | m**3 | kmol |
* | 3 | VCS_UNITS_MKS | Joules / Kmol (Cantera) | Pa | m**3 | kmol |
*
* see vcs_defs.h for more information.
*
* Currently, this value should be the same as the owning VCS_PROB or
* VCS_SOLVE object. There is no code for handling anything else atm.
*
* (This variable is needed for the vcsc code, where it is not equal to
* VCS_UNITS_MKS).
*/
int p_VCS_UnitsFormat;
//! Convention for the activity formulation
/*!
* * 0 = molar based activities (default)
@ -773,9 +745,7 @@ private:
/*!
* Note, This is the chemical potential derived strictly from the polynomial
* in temperature. Pressure effects have to be added in to get to the
* standard state.
*
* Units -> depends on VCS_UnitsFormat variable. Cantera -> J/kmol
* standard state. Units are J/kmol.
*/
mutable vector_fp SS0ChemicalPotential;
@ -783,9 +753,7 @@ private:
//! current Temperature and pressure.
/*!
* Note, This is the chemical potential at unit activity. Thus, we can call
* it the standard state chemical potential as well.
*
* Units -> depends on VCS_UnitsFormat variable. Cantera -> J/kmol.
* it the standard state chemical potential as well. Units are J/kmol.
*/
mutable vector_fp StarChemicalPotential;

35
include/cantera/equil/vcs_defs.h
View file

@ -227,41 +227,6 @@ namespace Cantera
//@}
/*!
* @name Units for the chemical potential data and pressure variables
*
* @verbatim
Chem_Pot Pres vol moles
-------------------------------------------------
VCS_UNITS_KCALMOL = kcal/mol Pa m**3 kmol
VCS_UNITS_UNITLESS = MU / RT -> no units Pa m**3 kmol
VCS_UNITS_KJMOL = kJ / mol Pa m**3 kmol
VCS_UNITS_KELVIN = KELVIN -> MU / R Pa m**3 kmol
VCS_UNITS_MKS = Joules / Kmol (Cantera) Pa m**3 kmol
Energy:
VCS_UNITS_KCALMOL = kcal/mol
VCS_UNITS_UNITLESS = MU / RT -> no units
VCS_UNITS_KJMOL = kJ / mol
VCS_UNITS_KELVIN = KELVIN -> MU / R
VCS_UNITS_MKS = J / kmol
Pressure: (Pref and Pres)
VCS_UNITS_KCALMOL = Pa
VCS_UNITS_UNITLESS = Pa
VCS_UNITS_KJMOL = Pa
VCS_UNITS_KELVIN = Pa
VCS_UNITS_MKS = Pa = kg / m s2
@endverbatim
* @{
*/
#define VCS_UNITS_KCALMOL -1
#define VCS_UNITS_UNITLESS 0
#define VCS_UNITS_KJMOL 1
#define VCS_UNITS_KELVIN 2
#define VCS_UNITS_MKS 3
//@}
/*!
* @name Types of Element Constraint Equations
*

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

@ -86,17 +86,6 @@ public:
double T_Time_vcs;
};
//! Returns the value of the gas constant in the units specified by parameter
/*!
* @param mu_units Specifies the units.
* - VCS_UNITS_KCALMOL: kcal gmol-1 K-1
* - VCS_UNITS_UNITLESS: 1.0 K-1
* - VCS_UNITS_KJMOL: kJ gmol-1 K-1
* - VCS_UNITS_KELVIN: 1.0 K-1
* - VCS_UNITS_MKS: joules kmol-1 K-1 = kg m2 s-2 kmol-1 K-1
*/
double vcsUtil_gasConstant(int mu_units);
//! Definition of the function pointer for the root finder
/*!
* see vcsUtil_root1d for a definition of how to use this.

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

@ -48,7 +48,7 @@ public:
size_t NPHASE0;
//! Vector of chemical potentials of the species. This is a calculated
//! output quantity. length = number of species. units = m_VCS_UnitsFormat
//! output quantity. length = number of species.
vector_fp m_gibbsSpecies;
//! Total number of moles of the kth species.
@ -59,8 +59,6 @@ public:
*
* On output, this contains the solution for the total number of moles of
* the kth species.
*
* units = m_VCS_UnitsFormat
*/
vector_fp w;
@ -71,7 +69,7 @@ public:
//! Element abundances for jth element
/*!
* This is input from the input file and is considered a constant from
* thereon within the vcs_solve_TP(). units = m_VCS_UnitsFormat
* thereon within the vcs_solve_TP().
*/
vector_fp gai;
@ -101,15 +99,10 @@ public:
double T;
//! Pressure
/*!
* units given by m_VCS_UnitsFormat
* -> are now PA
*/
double PresPA;
//! Volume of the entire system
/*!
* units given by m_VCS_UnitsFormat
* Note, this is an output variable atm
*/
double Vol;
@ -121,24 +114,6 @@ public:
*/
vector_fp VolPM;
//! Units for the chemical potential data, pressure data, volume, and
//! species amounts
/*!
* All internally stored quantities will have these units. Also, printed
* quantities will display in these units.
*
* | | | Chem_Pot | Pres | vol | moles|
* |---|----------------------|-------------------------|------|-------|------|
* |-1 | `VCS_UNITS_KCALMOL` | kcal/mol | atm | cm**3 | gmol |
* | 0 | `VCS_UNITS_UNITLESS` | MU / RT -> no units | atm | cm**3 | gmol |
* | 1 | `VCS_UNITS_KJMOL` | kJ / mol | atm | cm**3 | gmol |
* | 2 | `VCS_UNITS_KELVIN` | KELVIN -> MU / R | atm | cm**3 | gmol |
* | 3 | `VCS_UNITS_MKS` | Joules / Kmol (Cantera) | Pa | m**3 | kmol |
*
* see vcs_defs.h for more information
*/
int m_VCS_UnitsFormat;
//! Specification of the initial estimate method
/*!
* * 0: user estimate

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

@ -646,7 +646,7 @@ public:
* 0 -> Dont print any intermediate results
* @param maxit Maximum number of iterations for the algorithm
* @param T Value of the Temperature (Kelvin)
* @param pres Value of the Pressure (units given by m_VCS_UnitsFormat variable
* @param pres Value of the Pressure
* @return Returns an integer representing the success of the algorithm
* * 0 = Equilibrium Achieved
* * 1 = Range space error encountered. The element abundance criteria are
@ -893,26 +893,6 @@ public:
*/
int vcs_rearrange();
//! Returns the multiplier for electric charge terms
/*
* This is basically equal to F/RT
*
* @param mu_units integer representing the dimensional units system
* @param TKelvin double Temperature in Kelvin
* @returns the value of F/RT
*/
double vcs_nondim_Farad(int mu_units, double TKelvin) const;
//! Returns the multiplier for the nondimensionalization of the equations
/*!
* This is basically equal to RT
*
* @param mu_units integer representing the dimensional units system
* @param TKelvin double Temperature in Kelvin
* @returns the value of RT
*/
double vcs_nondimMult_TP(int mu_units, double TKelvin) const;
//! Nondimensionalize the problem data
/*!
* Nondimensionalize the free energies using the divisor, R * T
@ -938,14 +918,6 @@ public:
*/
void vcs_redim_TP();
//! Print the string representing the Chemical potential units
/*!
* This gets printed using plogf()
*
* @param unitsFormat Integer representing the units system
*/
void vcs_printChemPotUnits(int unitsFormat) const;
//! Computes the current elemental abundances vector
/*!
* Computes the elemental abundances vector, m_elemAbundances[], and stores
@ -1636,18 +1608,7 @@ public:
//! Temperature (Kelvin)
double m_temperature;
//! Pressure (units are determined by m_VCS_UnitsFormat
/*!
* | Values | units |
* | ------ | -----
* | -1: | atm |
* | 0: | atm |
* | 1: | atm |
* | 2: | atm |
* | 3: | Pa |
*
* Units being changed to Pa
*/
//! Pressure
double m_pressurePA;
//! Total kmoles of inert to add to each phase
@ -1869,7 +1830,7 @@ public:
//! Partial molar volumes of the species
/*!
* units = mks (m^3/kmol) -determined by m_VCS_UnitsFormat
* units = mks (m^3/kmol)
* Length = number of species
*/
vector_fp m_PMVolumeSpecies;
@ -1901,18 +1862,6 @@ public:
*/
int m_timing_print_lvl;
//! Units for the chemical potential data
/*!
* | Value | chemical potential units | pressure units |
* | ----- | ------------------------ | -------------- |
* | -1 | kcal/mol | Pa |
* | 0 | MU/RT | Pa |
* | 1 | kJ/mol | Pa |
* | 2 | Kelvin | Pa |
* | 3 | J / kmol | Pa |
*/
int m_VCS_UnitsFormat;
friend class vcs_phaseStabilitySolve;
};

10
include/cantera/equil/vcs_species_thermo.h
View file

@ -30,8 +30,7 @@ class vcs_VolPhase;
/*!
* Identifies the thermo model for the species. This structure is shared by
* volumetric and surface species. However, each will have its own types of
* thermodynamic models. These quantities all have appropriate units. The units
* are specified by VCS_UnitsFormat.
* thermodynamic models. These quantities all have appropriate units.
*/
class VCS_SPECIES_THERMO
{
@ -85,8 +84,6 @@ public:
//! parameter that is used in the VCS_SSVOL_CONSTANT model.
double SSStar_Vol0;
int m_VCS_UnitsFormat;
VCS_SPECIES_THERMO(size_t indexPhase, size_t indexSpeciesPhase);
virtual ~VCS_SPECIES_THERMO() {}
@ -102,7 +99,7 @@ public:
*
* @param kspec species global index
* @param TKelvin Temperature in Kelvin
* @param pres pressure is given in units specified by if__ variable.
* @param pres pressure in Pa
* @return standard state free energy in units of Kelvin.
*/
virtual double GStar_R_calc(size_t kspec, double TKelvin, double pres);
@ -121,8 +118,7 @@ public:
* This function calculates the standard state molar volume for species,
* kspec, at the temperature TKelvin and pressure, Pres,
*
* @return standard state volume in cm**3 per mol.
* (if__=3) m**3 / kmol
* @return standard state volume in m**3 / kmol
*/
virtual double VolStar_calc(size_t kglob, double TKelvin, double Pres);

46
src/equil/vcs_MultiPhaseEquil.cpp
View file

@ -514,24 +514,8 @@ int vcs_MultiPhaseEquil::equilibrate_TP(int estimateEquil,
plogf("\n");
plogf("----------------------------------------"
"---------------------\n");
plogf(" Name Mole_Number");
if (m_vprob.m_VCS_UnitsFormat == VCS_UNITS_MKS) {
plogf("(kmol)");
} else {
plogf("(gmol)");
}
plogf(" Mole_Fraction Chem_Potential");
if (m_vprob.m_VCS_UnitsFormat == VCS_UNITS_KCALMOL) {
plogf(" (kcal/mol)\n");
} else if (m_vprob.m_VCS_UnitsFormat == VCS_UNITS_UNITLESS) {
plogf(" (Dimensionless)\n");
} else if (m_vprob.m_VCS_UnitsFormat == VCS_UNITS_KJMOL) {
plogf(" (kJ/mol)\n");
} else if (m_vprob.m_VCS_UnitsFormat == VCS_UNITS_KELVIN) {
plogf(" (Kelvin)\n");
} else if (m_vprob.m_VCS_UnitsFormat == VCS_UNITS_MKS) {
plogf(" (J/kmol)\n");
}
plogf(" Name Mole_Number(kmol)");
plogf(" Mole_Fraction Chem_Potential (J/kmol)\n");
plogf("--------------------------------------------------"
"-----------\n");
for (size_t i = 0; i < m_vprob.nspecies; i++) {
@ -713,7 +697,6 @@ int vcs_Cantera_to_vprob(MultiPhase* mphase, VCS_PROB* vprob)
vprob->nspecies = totNumSpecies;
vprob->ne = 0;
vprob->NPhase = totNumPhases;
vprob->m_VCS_UnitsFormat = VCS_UNITS_MKS;
// Set the initial estimate to a machine generated estimate for now
// We will work out the details later.
vprob->iest = -1;
@ -752,7 +735,6 @@ int vcs_Cantera_to_vprob(MultiPhase* mphase, VCS_PROB* vprob)
VolPhase->m_gasPhase = gasPhase;
// Tell the vcs_VolPhase pointer about cantera
VolPhase->p_VCS_UnitsFormat = vprob->m_VCS_UnitsFormat;
VolPhase->setPtrThermoPhase(tPhase);
VolPhase->setTotalMoles(0.0);
@ -874,7 +856,6 @@ int vcs_Cantera_to_vprob(MultiPhase* mphase, VCS_PROB* vprob)
// Transfer the thermo specification of the species
// vprob->SpeciesThermo[]
ts_ptr->m_VCS_UnitsFormat = VolPhase->p_VCS_UnitsFormat;
// Add lookback connectivity into the thermo object first
ts_ptr->IndexPhase = iphase;
@ -945,11 +926,10 @@ int vcs_Cantera_to_vprob(MultiPhase* mphase, VCS_PROB* vprob)
// Now, calculate a sample naught Gibbs free energy calculation
// at the specified temperature.
double R = vcsUtil_gasConstant(vprob->m_VCS_UnitsFormat);
for (size_t k = 0; k < nSpPhase; k++) {
vcs_SpeciesProperties* sProp = VolPhase->speciesProperty(k);
ts_ptr = sProp->SpeciesThermo;
ts_ptr->SS0_feSave = VolPhase->G0_calc_one(k)/ R;
ts_ptr->SS0_feSave = VolPhase->G0_calc_one(k)/ GasConstant;
ts_ptr->SS0_TSave = vprob->T;
}
}
@ -1214,24 +1194,8 @@ int vcs_MultiPhaseEquil::determine_PhaseStability(int iph, double& funcStab, int
plogf("\n");
plogf("----------------------------------------"
"---------------------\n");
plogf(" Name Mole_Number");
if (m_vprob.m_VCS_UnitsFormat == VCS_UNITS_MKS) {
plogf("(kmol)");
} else {
plogf("(gmol)");
}
plogf(" Mole_Fraction Chem_Potential");
if (m_vprob.m_VCS_UnitsFormat == VCS_UNITS_KCALMOL) {
plogf(" (kcal/mol)\n");
} else if (m_vprob.m_VCS_UnitsFormat == VCS_UNITS_UNITLESS) {
plogf(" (Dimensionless)\n");
} else if (m_vprob.m_VCS_UnitsFormat == VCS_UNITS_KJMOL) {
plogf(" (kJ/mol)\n");
} else if (m_vprob.m_VCS_UnitsFormat == VCS_UNITS_KELVIN) {
plogf(" (Kelvin)\n");
} else if (m_vprob.m_VCS_UnitsFormat == VCS_UNITS_MKS) {
plogf(" (J/kmol)\n");
}
plogf(" Name Mole_Number(kmol)");
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]);

12
src/equil/vcs_TP.cpp
View file

@ -45,24 +45,12 @@ int VCS_SOLVE::vcs_TP(int ipr, int ip1, int maxit, double T_arg, double pres_arg
int VCS_SOLVE::vcs_evalSS_TP(int ipr, int ip1, double Temp, double pres)
{
// We need to special case VCS_UNITS_UNITLESS, here. cpc_ts_GStar_calc()
// returns units of Kelvin. Also, the temperature comes into play in
// calculating the ideal equation of state contributions, and other
// equations of state also. Therefore, we will emulate the VCS_UNITS_KELVIN
// case, here by changing the initial Gibbs free energy units to Kelvin
// before feeding them to the cpc_ts_GStar_calc() routine. Then, we will
// revert them back to unitless at the end of this routine.
for (size_t iph = 0; iph < m_numPhases; iph++) {
vcs_VolPhase* vph = m_VolPhaseList[iph];
vph->setState_TP(m_temperature, m_pressurePA);
vph->sendToVCS_GStar(&m_SSfeSpecies[0]);
}
if (m_VCS_UnitsFormat == VCS_UNITS_UNITLESS) {
for (size_t i = 0; i < m_numSpeciesTot; ++i) {
m_SSfeSpecies[i] /= Temp;
}
}
return VCS_SUCCESS;
}

4
src/equil/vcs_VolPhase.cpp
View file

@ -26,7 +26,6 @@ vcs_VolPhase::vcs_VolPhase(VCS_SOLVE* owningSolverObject) :
m_gasPhase(false),
m_eqnState(VCS_EOS_CONSTANT),
ChargeNeutralityElement(npos),
p_VCS_UnitsFormat(VCS_UNITS_MKS),
p_activityConvention(0),
m_numElemConstraints(0),
m_elemGlobalIndex(0),
@ -68,7 +67,6 @@ vcs_VolPhase::vcs_VolPhase(const vcs_VolPhase& b) :
m_gasPhase(b.m_gasPhase),
m_eqnState(b.m_eqnState),
ChargeNeutralityElement(b.ChargeNeutralityElement),
p_VCS_UnitsFormat(b.p_VCS_UnitsFormat),
p_activityConvention(b.p_activityConvention),
m_numElemConstraints(b.m_numElemConstraints),
m_numSpecies(b.m_numSpecies),
@ -113,7 +111,6 @@ vcs_VolPhase& vcs_VolPhase::operator=(const vcs_VolPhase& b)
m_gasPhase = b.m_gasPhase;
m_eqnState = b.m_eqnState;
ChargeNeutralityElement = b.ChargeNeutralityElement;
p_VCS_UnitsFormat = b.p_VCS_UnitsFormat;
p_activityConvention= b.p_activityConvention;
m_numSpecies = b.m_numSpecies;
m_numElemConstraints = b.m_numElemConstraints;
@ -715,7 +712,6 @@ void vcs_VolPhase::setPtrThermoPhase(ThermoPhase* tp_ptr)
Temp_ = TP_ptr->temperature();
Pres_ = TP_ptr->pressure();
setState_TP(Temp_, Pres_);
p_VCS_UnitsFormat = VCS_UNITS_MKS;
m_phi = TP_ptr->electricPotential();
size_t nsp = TP_ptr->nSpecies();
size_t nelem = TP_ptr->nElements();

142
src/equil/vcs_nondim.cpp
View file

@ -15,51 +15,12 @@
namespace Cantera
{
double VCS_SOLVE::vcs_nondim_Farad(int mu_units, double TKelvin) const
{
if (TKelvin <= 0.0) {
TKelvin = 293.15;
}
switch (mu_units) {
case VCS_UNITS_MKS:
case VCS_UNITS_KJMOL:
case VCS_UNITS_KCALMOL:
return ElectronCharge * Avogadro / (TKelvin * GasConstant);
case VCS_UNITS_UNITLESS:
return ElectronCharge * Avogadro;
case VCS_UNITS_KELVIN:
return ElectronCharge * Avogadro/ TKelvin;
default:
throw CanteraError("vcs_nondim_Farad", "unknown units: {}", mu_units);
}
}
double VCS_SOLVE::vcs_nondimMult_TP(int mu_units, double TKelvin) const
{
if (TKelvin <= 0.0) {
TKelvin = 293.15;
}
switch (mu_units) {
case VCS_UNITS_KCALMOL:
return TKelvin * GasConst_cal_mol_K * 1e-3;
case VCS_UNITS_UNITLESS:
return 1.0;
case VCS_UNITS_KJMOL:
return TKelvin * GasConstant * 1e-6;
case VCS_UNITS_KELVIN:
return TKelvin;
case VCS_UNITS_MKS:
return TKelvin * GasConstant;
default:
throw CanteraError("vcs_nondimMult_TP", "unknown units: {}", mu_units);
}
}
void VCS_SOLVE::vcs_nondim_TP()
{
if (m_unitsState == VCS_DIMENSIONAL_G) {
m_unitsState = VCS_NONDIMENSIONAL_G;
double tf = 1.0 / vcs_nondimMult_TP(m_VCS_UnitsFormat, m_temperature);
double tf = 1.0 / (GasConstant * m_temperature);
for (size_t i = 0; i < m_numSpeciesTot; ++i) {
// Modify the standard state and total chemical potential data,
// FF(I), to make it dimensionless, i.e., mu / RT. Thus, we may
@ -70,7 +31,7 @@ void VCS_SOLVE::vcs_nondim_TP()
m_feSpecies_old[i] *= tf;
}
m_Faraday_dim = vcs_nondim_Farad(m_VCS_UnitsFormat, m_temperature);
m_Faraday_dim = ElectronCharge * Avogadro / (m_temperature * GasConstant);
// Scale the total moles if necessary: First find out the total moles
double tmole_orig = vcs_tmoles();
@ -104,26 +65,24 @@ void VCS_SOLVE::vcs_nondim_TP()
}
if (m_totalMoleScale != 1.0) {
if (m_VCS_UnitsFormat == VCS_UNITS_MKS) {
if (m_debug_print_lvl >= 2) {
plogf(" --- vcs_nondim_TP() called: USING A MOLE SCALE OF %g until further notice", m_totalMoleScale);
plogendl();
}
for (size_t i = 0; i < m_numSpeciesTot; ++i) {
if (m_speciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
m_molNumSpecies_old[i] *= (1.0 / m_totalMoleScale);
}
}
for (size_t i = 0; i < m_numElemConstraints; ++i) {
m_elemAbundancesGoal[i] *= (1.0 / m_totalMoleScale);
if (m_debug_print_lvl >= 2) {
plogf(" --- vcs_nondim_TP() called: USING A MOLE SCALE OF %g until further notice", m_totalMoleScale);
plogendl();
}
for (size_t i = 0; i < m_numSpeciesTot; ++i) {
if (m_speciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
m_molNumSpecies_old[i] *= (1.0 / m_totalMoleScale);
}
}
for (size_t i = 0; i < m_numElemConstraints; ++i) {
m_elemAbundancesGoal[i] *= (1.0 / m_totalMoleScale);
}
for (size_t iph = 0; iph < m_numPhases; iph++) {
TPhInertMoles[iph] *= (1.0 / m_totalMoleScale);
if (TPhInertMoles[iph] != 0.0) {
vcs_VolPhase* vphase = m_VolPhaseList[iph];
vphase->setTotalMolesInert(TPhInertMoles[iph]);
}
for (size_t iph = 0; iph < m_numPhases; iph++) {
TPhInertMoles[iph] *= (1.0 / m_totalMoleScale);
if (TPhInertMoles[iph] != 0.0) {
vcs_VolPhase* vphase = m_VolPhaseList[iph];
vphase->setTotalMolesInert(TPhInertMoles[iph]);
}
}
vcs_tmoles();
@ -135,7 +94,7 @@ void VCS_SOLVE::vcs_redim_TP()
{
if (m_unitsState != VCS_DIMENSIONAL_G) {
m_unitsState = VCS_DIMENSIONAL_G;
double tf = vcs_nondimMult_TP(m_VCS_UnitsFormat, m_temperature);
double tf = m_temperature * GasConstant;
for (size_t i = 0; i < m_numSpeciesTot; ++i) {
// Modify the standard state and total chemical potential data,
@ -148,52 +107,27 @@ void VCS_SOLVE::vcs_redim_TP()
m_Faraday_dim *= tf;
}
if (m_totalMoleScale != 1.0) {
if (m_VCS_UnitsFormat == VCS_UNITS_MKS) {
if (m_debug_print_lvl >= 2) {
plogf(" --- vcs_redim_TP() called: getting rid of mole scale of %g", m_totalMoleScale);
plogendl();
}
for (size_t i = 0; i < m_numSpeciesTot; ++i) {
if (m_speciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
m_molNumSpecies_old[i] *= m_totalMoleScale;
}
}
for (size_t i = 0; i < m_numElemConstraints; ++i) {
m_elemAbundancesGoal[i] *= m_totalMoleScale;
}
for (size_t iph = 0; iph < m_numPhases; iph++) {
TPhInertMoles[iph] *= m_totalMoleScale;
if (TPhInertMoles[iph] != 0.0) {
vcs_VolPhase* vphase = m_VolPhaseList[iph];
vphase->setTotalMolesInert(TPhInertMoles[iph]);
}
}
vcs_tmoles();
if (m_debug_print_lvl >= 2) {
plogf(" --- vcs_redim_TP() called: getting rid of mole scale of %g", m_totalMoleScale);
plogendl();
}
for (size_t i = 0; i < m_numSpeciesTot; ++i) {
if (m_speciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
m_molNumSpecies_old[i] *= m_totalMoleScale;
}
}
for (size_t i = 0; i < m_numElemConstraints; ++i) {
m_elemAbundancesGoal[i] *= m_totalMoleScale;
}
}
}
void VCS_SOLVE::vcs_printChemPotUnits(int unitsFormat) const
{
switch (unitsFormat) {
case VCS_UNITS_KCALMOL:
plogf("kcal/gmol");
break;
case VCS_UNITS_UNITLESS:
plogf("dimensionless");
break;
case VCS_UNITS_KJMOL:
plogf("kJ/gmol");
break;
case VCS_UNITS_KELVIN:
plogf("Kelvin");
break;
case VCS_UNITS_MKS:
plogf("J/kmol");
break;
default:
throw CanteraError("VCS_SOLVE::vcs_printChemPotUnits", "unknown units!");
for (size_t iph = 0; iph < m_numPhases; iph++) {
TPhInertMoles[iph] *= m_totalMoleScale;
if (TPhInertMoles[iph] != 0.0) {
vcs_VolPhase* vphase = m_VolPhaseList[iph];
vphase->setTotalMolesInert(TPhInertMoles[iph]);
}
}
vcs_tmoles();
}
}

22
src/equil/vcs_prob.cpp
View file

@ -34,8 +34,6 @@ VCS_PROB::VCS_PROB(size_t nsp, size_t nel, size_t nph) :
T(298.15),
PresPA(1.0),
Vol(0.0),
// Set the units for the chemical potential data to be unitless
m_VCS_UnitsFormat(VCS_UNITS_UNITLESS),
// The default is to not expect an initial estimate of the species
// concentrations
iest(-1),
@ -224,30 +222,14 @@ void VCS_PROB::prob_report(int print_lvl)
plogf("\nElemental Abundances: ");
plogf(" Target_kmol ElemType ElActive\n");
double fac = 1.0;
if (m_VCS_UnitsFormat == VCS_UNITS_MKS) {
fac = 1.0;
}
for (size_t i = 0; i < ne; ++i) {
writeline(' ', 26, false);
plogf("%-2.2s", ElName[i]);
plogf("%20.12E ", fac * gai[i]);
plogf("%20.12E ", gai[i]);
plogf("%3d %3d\n", m_elType[i], ElActive[i]);
}
plogf("\nChemical Potentials: ");
if (m_VCS_UnitsFormat == VCS_UNITS_UNITLESS) {
plogf("(unitless)");
} else if (m_VCS_UnitsFormat == VCS_UNITS_KCALMOL) {
plogf("(kcal/gmol)");
} else if (m_VCS_UnitsFormat == VCS_UNITS_KJMOL) {
plogf("(kJ/gmol)");
} else if (m_VCS_UnitsFormat == VCS_UNITS_KELVIN) {
plogf("(Kelvin)");
} else if (m_VCS_UnitsFormat == VCS_UNITS_MKS) {
plogf("(J/kmol)");
}
plogf("\n");
plogf("\nChemical Potentials: (J/kmol)\n");
plogf(" Species (phase) "
" SS0ChemPot StarChemPot\n");
for (size_t iphase = 0; iphase < NPhase; iphase++) {

5
src/equil/vcs_report.cpp
View file

@ -252,10 +252,7 @@ int VCS_SOLVE::vcs_report(int iconv)
writeline('-', 93, true, true);
plogf("Chemical Potentials of the Species: (dimensionless)\n");
double rt = vcs_nondimMult_TP(m_VCS_UnitsFormat, m_temperature);
plogf("\t\t(RT = %g ", rt);
vcs_printChemPotUnits(m_VCS_UnitsFormat);
plogf(")\n");
plogf("\t\t(RT = %g J/kmol)\n", GasConstant * m_temperature);
plogf(" Name TKMoles StandStateChemPot "
" ln(AC) ln(X_i) | F z_i phi | ChemPot | (-lnMnaught)");
plogf("| (MolNum ChemPot)|");

9
src/equil/vcs_solve.cpp
View file

@ -49,8 +49,7 @@ VCS_SOLVE::VCS_SOLVE() :
m_Faraday_dim(ElectronCharge * Avogadro),
m_VCount(0),
m_debug_print_lvl(0),
m_timing_print_lvl(1),
m_VCS_UnitsFormat(VCS_UNITS_UNITLESS)
m_timing_print_lvl(1)
{
}
@ -88,8 +87,6 @@ void VCS_SOLVE::vcs_initSizes(const size_t nspecies0, const size_t nelements,
" Number of species is nonpositive\n");
}
m_VCS_UnitsFormat = VCS_UNITS_UNITLESS;
/*
* We will initialize sc[] to note the fact that it needs to be
* filled with meaningful information.
@ -464,9 +461,6 @@ int VCS_SOLVE::vcs_prob_specifyFully(const VCS_PROB* pub)
TPhInertMoles[iph] = Vphase->totalMolesInert();
}
// if__ : Copy over the units for the chemical potential
m_VCS_UnitsFormat = pub->m_VCS_UnitsFormat;
// tolerance requirements -> copy them over here and later
m_tolmaj = pub->tolmaj;
m_tolmin = pub->tolmin;
@ -629,7 +623,6 @@ int VCS_SOLVE::vcs_prob_specify(const VCS_PROB* pub)
m_temperature = pub->T;
m_pressurePA = pub->PresPA;
m_VCS_UnitsFormat = pub->m_VCS_UnitsFormat;
m_doEstimateEquil = pub->iest;
m_totalVol = pub->Vol;
m_tolmaj = pub->tolmaj;

18
src/equil/vcs_solve_TP.cpp
View file

@ -111,21 +111,7 @@ int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit)
if (m_doEstimateEquil == 0) {
plogf("\n USER ESTIMATE OF EQUILIBRIUM\n");
}
if (m_VCS_UnitsFormat == VCS_UNITS_KCALMOL) {
plogf(" Stan. Chem. Pot. in kcal/mole\n");
}
if (m_VCS_UnitsFormat == VCS_UNITS_UNITLESS) {
plogf(" Stan. Chem. Pot. is MU/RT\n");
}
if (m_VCS_UnitsFormat == VCS_UNITS_KJMOL) {
plogf(" Stan. Chem. Pot. in KJ/mole\n");
}
if (m_VCS_UnitsFormat == VCS_UNITS_KELVIN) {
plogf(" Stan. Chem. Pot. in Kelvin\n");
}
if (m_VCS_UnitsFormat == VCS_UNITS_MKS) {
plogf(" Stan. Chem. Pot. in J/kmol\n");
}
plogf(" Stan. Chem. Pot. in J/kmol\n");
plogf("\n SPECIES FORMULA VECTOR ");
writeline(' ', 41, false);
plogf(" STAN_CHEM_POT EQUILIBRIUM_EST. Species_Type\n\n");
@ -134,7 +120,7 @@ int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit)
plogf("%-4.4s ", m_elementName[i]);
}
plogf(" PhaseID\n");
double RT = vcs_nondimMult_TP(m_VCS_UnitsFormat, m_temperature);
double RT = GasConstant * m_temperature;
for (size_t i = 0; i < m_numSpeciesTot; ++i) {
plogf(" %-18.18s", m_speciesName[i]);
for (size_t j = 0; j < m_numElemConstraints; ++j) {

25
src/equil/vcs_species_thermo.cpp
View file

@ -33,8 +33,7 @@ VCS_SPECIES_THERMO::VCS_SPECIES_THERMO(size_t indexPhase,
SS0_Pref(1.01325E5),
SSStar_Model(VCS_SSSTAR_CONSTANT),
SSStar_Vol_Model(VCS_SSVOL_IDEALGAS),
SSStar_Vol0(-1.0),
m_VCS_UnitsFormat(VCS_UNITS_UNITLESS)
SSStar_Vol0(-1.0)
{
SS0_Pref = 1.01325E5;
}
@ -53,8 +52,7 @@ VCS_SPECIES_THERMO::VCS_SPECIES_THERMO(const VCS_SPECIES_THERMO& b) :
SS0_Pref(b.SS0_Pref),
SSStar_Model(b.SSStar_Model),
SSStar_Vol_Model(b.SSStar_Vol_Model),
SSStar_Vol0(b.SSStar_Vol0),
m_VCS_UnitsFormat(b.m_VCS_UnitsFormat)
SSStar_Vol0(b.SSStar_Vol0)
{
}
@ -76,7 +74,6 @@ VCS_SPECIES_THERMO::operator=(const VCS_SPECIES_THERMO& b)
SSStar_Model = b.SSStar_Model;
SSStar_Vol_Model = b.SSStar_Vol_Model;
SSStar_Vol0 = b.SSStar_Vol0;
m_VCS_UnitsFormat = b.m_VCS_UnitsFormat;
}
return *this;
}
@ -90,23 +87,14 @@ double VCS_SPECIES_THERMO::GStar_R_calc(size_t kglob, double TKelvin,
double pres)
{
double fe = G0_R_calc(kglob, TKelvin);
if (m_VCS_UnitsFormat != VCS_UNITS_MKS) {
throw CanteraError("VCS_SPECIES_THERMO::GStar_R_calc",
"Possible inconsistency");
}
OwningPhase->setState_TP(TKelvin, pres);
fe = OwningPhase->GStar_calc_one(IndexSpeciesPhase);
double R = vcsUtil_gasConstant(m_VCS_UnitsFormat);
return fe / R;
return fe / GasConstant;
}
double VCS_SPECIES_THERMO::VolStar_calc(size_t kglob, double TKelvin,
double presPA)
{
if (m_VCS_UnitsFormat != VCS_UNITS_MKS) {
throw CanteraError("VCS_SPECIES_THERMO::VolStar_calc",
"Possible inconsistency");
}
OwningPhase->setState_TP(TKelvin, presPA);
return OwningPhase->VolStar_calc_one(IndexSpeciesPhase);
}
@ -119,14 +107,9 @@ double VCS_SPECIES_THERMO::G0_R_calc(size_t kglob, double TKelvin)
if (TKelvin == SS0_TSave) {
return SS0_feSave;
}
if (m_VCS_UnitsFormat != VCS_UNITS_MKS) {
throw CanteraError("VCS_SPECIES_THERMO::G0_R_calc",
"Possible inconsistency");
}
OwningPhase->setState_T(TKelvin);
double fe = OwningPhase->G0_calc_one(IndexSpeciesPhase);
double R = vcsUtil_gasConstant(m_VCS_UnitsFormat);
fe /= R;
fe /= GasConstant;
SS0_feSave = fe;
SS0_TSave = TKelvin;
return fe;

19
src/equil/vcs_util.cpp
View file

@ -57,25 +57,6 @@ size_t vcs_optMax(const double* x, const double* xSize, size_t j, size_t n)
return largest;
}
double vcsUtil_gasConstant(int mu_units)
{
switch (mu_units) {
case VCS_UNITS_KCALMOL:
return GasConst_cal_mol_K * 1e-3;
case VCS_UNITS_UNITLESS:
return 1.0;
case VCS_UNITS_KJMOL:
return GasConstant * 1e-6;
case VCS_UNITS_KELVIN:
return 1.0;
case VCS_UNITS_MKS:
// joules / kg-mol K = kg m2 / s2 kg-mol K
return GasConstant;
default:
throw CanteraError("vcsUtil_gasConstant", "uknown units: {}", mu_units);
}
}
const char* vcs_speciesType_string(int speciesStatus, int length)
{
switch (speciesStatus) {