Internal rearrangement of vcs_VolPhase to handle setState StateCalc concepts.

(incremental)
This commit is contained in:
Harry Moffat 2008-06-20 16:17:46 +00:00
parent 59481c2808
commit 2057b2007a
6 changed files with 172 additions and 184 deletions

View file

@ -1168,16 +1168,16 @@ namespace VCSnonideal {
*/
switch (eos) {
case cIdealGas:
VolPhase->EqnState = VCS_EOS_IDEAL_GAS;
VolPhase->m_eqnState = VCS_EOS_IDEAL_GAS;
break;
case cIncompressible:
VolPhase->EqnState = VCS_EOS_CONSTANT;
VolPhase->m_eqnState = VCS_EOS_CONSTANT;
break;
case cSurf:
plogf("cSurf not handled yet\n");
exit(-1);
case cStoichSubstance:
VolPhase->EqnState = VCS_EOS_STOICH_SUB;
VolPhase->m_eqnState = VCS_EOS_STOICH_SUB;
break;
case cPureFluid:
if (printLvl > 1) {
@ -1190,13 +1190,13 @@ namespace VCSnonideal {
case cIdealSolidSolnPhase0:
case cIdealSolidSolnPhase1:
case cIdealSolidSolnPhase2:
VolPhase->EqnState = VCS_EOS_IDEAL_SOLN;
VolPhase->m_eqnState = VCS_EOS_IDEAL_SOLN;
break;
default:
if (printLvl > 1) {
plogf("Unknown Cantera EOS to VCSnonideal: %d\n", eos);
}
VolPhase->EqnState = VCS_EOS_UNK_CANTERA;
VolPhase->m_eqnState = VCS_EOS_UNK_CANTERA;
if (!VolPhase->usingCanteraCalls()) {
plogf("vcs functions asked for, but unimplemented\n");
exit(-1);
@ -1440,7 +1440,7 @@ namespace VCSnonideal {
for (int iphase = 0; iphase < vprob->NPhase; iphase++) {
vcs_VolPhase *VolPhase = vprob->VPhaseList[iphase];
std::string sEOS = string16_EOSType(VolPhase->EqnState);
std::string sEOS = string16_EOSType(VolPhase->m_eqnState);
plogf("%16s %5d %5d %8d %16s %8d %16e ", VolPhase->PhaseName.c_str(),
VolPhase->VP_ID, VolPhase->SingleSpecies,
VolPhase->m_gasPhase, sEOS.c_str(),
@ -1558,7 +1558,7 @@ namespace VCSnonideal {
for (int iphase = 0; iphase < vprob->NPhase; iphase++) {
vcs_VolPhase *VolPhase = vprob->VPhaseList[iphase];
std::string sEOS = string16_EOSType(VolPhase->EqnState);
std::string sEOS = string16_EOSType(VolPhase->m_eqnState);
plogf("%16s %5d %5d %8d %16s %8d %16e ", VolPhase->PhaseName.c_str(),
VolPhase->VP_ID, VolPhase->SingleSpecies,
VolPhase->m_gasPhase, sEOS.c_str(),

View file

@ -35,7 +35,7 @@ namespace VCSnonideal {
Domain_ID(-1),
SingleSpecies(true),
m_gasPhase(false),
EqnState(VCS_EOS_CONSTANT),
m_eqnState(VCS_EOS_CONSTANT),
nElemConstraints(0),
ChargeNeutralityElement(-1),
ElGlobalIndex(0),
@ -99,7 +99,7 @@ namespace VCSnonideal {
Domain_ID(b.Domain_ID),
SingleSpecies(b.SingleSpecies),
m_gasPhase(b.m_gasPhase),
EqnState(b.EqnState),
m_eqnState(b.m_eqnState),
nElemConstraints(b.nElemConstraints),
ChargeNeutralityElement(b.ChargeNeutralityElement),
NVolSpecies(b.NVolSpecies),
@ -133,15 +133,14 @@ namespace VCSnonideal {
*/
*this = b;
}
/*****************************************************************************/
/***********************************************************************************/
/*
* Assignment operator()
*
* (note, this is used, so keep it current!)
*/
vcs_VolPhase& vcs_VolPhase::operator=(const vcs_VolPhase& b)
{
vcs_VolPhase& vcs_VolPhase::operator=(const vcs_VolPhase& b) {
int k;
if (&b != this) {
int old_num = NVolSpecies;
@ -155,7 +154,7 @@ namespace VCSnonideal {
Domain_ID = b.Domain_ID;
SingleSpecies = b.SingleSpecies;
m_gasPhase = b.m_gasPhase;
EqnState = b.EqnState;
m_eqnState = b.m_eqnState;
NVolSpecies = b.NVolSpecies;
nElemConstraints = b.nElemConstraints;
@ -334,15 +333,14 @@ namespace VCSnonideal {
m_UpToDate_VolPM = false;
m_UpToDate_GStar = false;
}
/*******************************************************************************/
/************************************************************************************/
//! Evaluate activity coefficients
/*!
* We carry out a calculation whenever UpTODate_AC is false. Specifically
* whenever a phase goes zero, we do not carry out calculations on it.
*/
void vcs_VolPhase::evaluateActCoeff() const {
if (m_UpToDate_AC == true) return;
void vcs_VolPhase::_updateActCoeff() const {
if (m_isIdealSoln) {
m_UpToDate_AC = true;
return;
@ -364,7 +362,7 @@ namespace VCSnonideal {
}
m_UpToDate_AC = true;
}
/********************************************************************************/
/***********************************************************************************/
/*
*
@ -374,7 +372,9 @@ namespace VCSnonideal {
* one.
*/
double vcs_VolPhase::AC_calc_one(int kspec) const {
evaluateActCoeff();
if (! m_UpToDate_AC) {
_updateActCoeff();
}
return(ActCoeff[kspec]);
}
/************************************************************************************/
@ -408,7 +408,7 @@ namespace VCSnonideal {
}
}
}
/***********************************************************************/
/*******************************************************************************/
// Gibbs free energy calculation at a temperature for the reference state
// of a species, return a value for one species
@ -422,7 +422,7 @@ namespace VCSnonideal {
G0_calc(tkelvin);
return SS0ChemicalPotential[kspec];
}
/***********************************************************************/
/*******************************************************************************/
// Gibbs free energy calculation for standard states
/*
@ -432,24 +432,22 @@ namespace VCSnonideal {
* @param TKelvin Current temperature
* @param pres Current pressure (pascal)
*/
void vcs_VolPhase::GStar_calc() const {
if (!m_UpToDate_GStar) {
if (m_useCanteraCalls) {
TP_ptr->getStandardChemPotentials(VCS_DATA_PTR(StarChemicalPotential));
} else {
double R = vcsUtil_gasConstant(m_VCS_UnitsFormat);
for (int k = 0; k < NVolSpecies; k++) {
int kglob = IndSpecies[k];
vcs_SpeciesProperties *sProp = ListSpeciesPtr[k];
VCS_SPECIES_THERMO *sTherm = sProp->SpeciesThermo;
StarChemicalPotential[k] =
R * (sTherm->GStar_R_calc(kglob, Temp, Pres));
}
void vcs_VolPhase::_updateGStar() const {
if (m_useCanteraCalls) {
TP_ptr->getStandardChemPotentials(VCS_DATA_PTR(StarChemicalPotential));
} else {
double R = vcsUtil_gasConstant(m_VCS_UnitsFormat);
for (int k = 0; k < NVolSpecies; k++) {
int kglob = IndSpecies[k];
vcs_SpeciesProperties *sProp = ListSpeciesPtr[k];
VCS_SPECIES_THERMO *sTherm = sProp->SpeciesThermo;
StarChemicalPotential[k] =
R * (sTherm->GStar_R_calc(kglob, Temp, Pres));
}
m_UpToDate_GStar = true;
}
m_UpToDate_GStar = true;
}
/***********************************************************************/
/*****************************************************************************/
// Gibbs free energy calculation for standard state of one species
/*
@ -465,11 +463,11 @@ namespace VCSnonideal {
*/
double vcs_VolPhase::GStar_calc_one(int kspec) {
if (!m_UpToDate_GStar) {
GStar_calc();
_updateGStar();
}
return StarChemicalPotential[kspec];
}
/***********************************************************************/
/*****************************************************************************/
// Set the mole fractions from a conventional mole fraction vector
/*
@ -540,7 +538,7 @@ namespace VCSnonideal {
if (molesSpeciesVCS == 0) {
#ifdef DEBUG_MODE
if (m_owningSolverObject == 0) {
printf("shouldn't be here\n");
printf("vcs_VolPhase::setMolesFromVCS shouldn't be here\n");
std::exit(-1);
}
#endif
@ -551,7 +549,7 @@ namespace VCSnonideal {
}
#ifdef DEBUG_MODE
else {
printf("shouldn't be here\n");
printf("vcs_VolPhase::setMolesFromVCS shouldn't be here\n");
std::exit(-1);
}
#endif
@ -561,12 +559,12 @@ namespace VCSnonideal {
if (m_owningSolverObject) {
if (stateCalc == VCS_STATECALC_OLD) {
if (molesSpeciesVCS != VCS_DATA_PTR(m_owningSolverObject->m_molNumSpecies_old)) {
printf("shouldn't be here\n");
printf("vcs_VolPhase::setMolesFromVCS shouldn't be here\n");
std::exit(-1);
}
} else if (stateCalc == VCS_STATECALC_NEW) {
if (molesSpeciesVCS != VCS_DATA_PTR(m_owningSolverObject->m_molNumSpecies_new)) {
printf("shouldn't be here\n");
printf("vcs_VolPhase::setMolesFromVCS shouldn't be here\n");
std::exit(-1);
}
}
@ -577,14 +575,16 @@ namespace VCSnonideal {
for (int k = 0; k < NVolSpecies; k++) {
if (SpeciesUnknownType[k] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
kglob = IndSpecies[k];
tmp = MAX(0.0, molesSpeciesVCS[kglob]);
Xmol[k] = tmp;
v_totalMoles += tmp;
v_totalMoles += MAX(0.0, molesSpeciesVCS[kglob]);
}
}
if (v_totalMoles > 0.0) {
for (int k = 0; k < NVolSpecies; k++) {
Xmol[k] /= v_totalMoles;
if (SpeciesUnknownType[k] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
kglob = IndSpecies[k];
tmp = MAX(0.0, molesSpeciesVCS[kglob]);
Xmol[k] = tmp / v_totalMoles;
}
}
Existence = 1;
} else {
@ -624,7 +624,7 @@ namespace VCSnonideal {
m_vcsStateStatus = stateCalc;
}
/***********************************************************************/
/**************************************************************************/
// Set the moles within the phase
/*
@ -657,7 +657,7 @@ namespace VCSnonideal {
}
}
}
/***********************************************************************/
/**************************************************************************/
// Update the moles within the phase, if necessary
/*
@ -680,7 +680,7 @@ namespace VCSnonideal {
}
}
}
/***********************************************************************/
/**************************************************************************/
// Fill in an activity coefficients vector within a VCS_SOLVE object
/*
@ -696,7 +696,7 @@ namespace VCSnonideal {
double * const AC) {
updateFromVCS_MoleNumbers(stateCalc);
if (!m_UpToDate_AC) {
evaluateActCoeff();
_updateActCoeff();
}
int kglob;
for (int k = 0; k < NVolSpecies; k++) {
@ -704,7 +704,7 @@ namespace VCSnonideal {
AC[kglob] = ActCoeff[k];
}
}
/***********************************************************************/
/****************************************************************************/
// Fill in the partial molar volume vector for VCS
/*
@ -718,7 +718,7 @@ namespace VCSnonideal {
*/
double vcs_VolPhase::sendToVCS_VolPM(double * const VolPM) const {
if (!m_UpToDate_VolPM) {
(void) VolPM_calc();
(void) _updateVolPM();
}
int kglob;
for (int k = 0; k < NVolSpecies; k++) {
@ -727,7 +727,7 @@ namespace VCSnonideal {
}
return m_totalVol;
}
/***********************************************************************/
/****************************************************************************/
// Fill in the partial molar volume vector for VCS
/*
@ -741,7 +741,7 @@ namespace VCSnonideal {
*/
void vcs_VolPhase::sendToVCS_GStar(double * const gstar){
if (!m_UpToDate_GStar) {
GStar_calc();
_updateGStar();
}
int kglob;
for (int k = 0; k < NVolSpecies; k++) {
@ -749,7 +749,7 @@ namespace VCSnonideal {
gstar[kglob] = StarChemicalPotential[k];
}
}
/***********************************************************************/
/****************************************************************************/
void vcs_VolPhase::setElectricPotential(double phi) {
@ -763,12 +763,12 @@ namespace VCSnonideal {
m_UpToDate_VolPM = false;
m_UpToDate_GStar = false;
}
/***********************************************************************/
/*****************************************************************************/
double vcs_VolPhase::electricPotential() const {
return m_phi;
}
/***********************************************************************/
/****************************************************************************/
// Sets the temperature and pressure in this object and
// underlying objects
@ -798,7 +798,7 @@ namespace VCSnonideal {
m_UpToDate_VolPM = false;
m_UpToDate_GStar = false;
}
/***********************************************************************/
/**************************************************************************/
// Molar volume calculation for standard states
/*
@ -808,26 +808,22 @@ namespace VCSnonideal {
* @param TKelvin Current temperature
* @param pres Current pressure (pascal)
*
* Calculations are in m**3/kmol
* Calculations are in m**3 / kmol
*/
void vcs_VolPhase::VolStar_calc() const {
if (!m_UpToDate_VolStar) {
if (m_useCanteraCalls) {
TP_ptr->getStandardVolumes(VCS_DATA_PTR(StarMolarVol));
} else {
for (int k = 0; k < NVolSpecies; k++) {
int kglob = IndSpecies[k];
vcs_SpeciesProperties *sProp = ListSpeciesPtr[k];
VCS_SPECIES_THERMO *sTherm = sProp->SpeciesThermo;
StarMolarVol[k] =
(sTherm->VolStar_calc(kglob, Temp, Pres));
}
void vcs_VolPhase::_updateVolStar() const {
if (m_useCanteraCalls) {
TP_ptr->getStandardVolumes(VCS_DATA_PTR(StarMolarVol));
} else {
for (int k = 0; k < NVolSpecies; k++) {
int kglob = IndSpecies[k];
vcs_SpeciesProperties *sProp = ListSpeciesPtr[k];
VCS_SPECIES_THERMO *sTherm = sProp->SpeciesThermo;
StarMolarVol[k] = (sTherm->VolStar_calc(kglob, Temp, Pres));
}
m_UpToDate_VolStar = true;
}
m_UpToDate_VolStar = true;
}
/***********************************************************************/
/*****************************************************************************/
// Molar volume calculation for standard state of one species
/*
@ -842,51 +838,51 @@ namespace VCSnonideal {
* @return molar volume of the kspec species's standard
* state
*/
double vcs_VolPhase::VolStar_calc_one(int kspec, double tkelvin,
double pres) {
setState_TP(tkelvin, pres);
double vcs_VolPhase::VolStar_calc_one(int kspec) const {
if (!m_UpToDate_VolStar) {
VolStar_calc();
_updateVolStar();
}
return StarMolarVol[kspec];
}
/****************************************************************************/
// Calculate the partial molar volumes of all species and return the
// total volume
/*
* Calculates these quantitites internally and then stores them
*
* VolPM_calc
* @return total volume (m**3)
*/
double vcs_VolPhase::VolPM_calc() const {
double vcs_VolPhase::_updateVolPM() const {
int k, kglob;
if (!m_UpToDate_VolPM) {
if (m_useCanteraCalls) {
TP_ptr->getPartialMolarVolumes(VCS_DATA_PTR(PartialMolarVol));
} else {
for (k = 0; k < NVolSpecies; k++) {
kglob = IndSpecies[k];
vcs_SpeciesProperties *sProp = ListSpeciesPtr[k];
VCS_SPECIES_THERMO *sTherm = sProp->SpeciesThermo;
StarMolarVol[k] = (sTherm->VolStar_calc(kglob, Temp, Pres));
}
for (k = 0; k < NVolSpecies; k++) {
PartialMolarVol[k] = StarMolarVol[k];
}
}
m_totalVol = 0.0;
if (m_useCanteraCalls) {
TP_ptr->getPartialMolarVolumes(VCS_DATA_PTR(PartialMolarVol));
} else {
for (k = 0; k < NVolSpecies; k++) {
m_totalVol += PartialMolarVol[k] * Xmol[k];
kglob = IndSpecies[k];
vcs_SpeciesProperties *sProp = ListSpeciesPtr[k];
VCS_SPECIES_THERMO *sTherm = sProp->SpeciesThermo;
StarMolarVol[k] = (sTherm->VolStar_calc(kglob, Temp, Pres));
}
m_totalVol *= v_totalMoles;
for (k = 0; k < NVolSpecies; k++) {
PartialMolarVol[k] = StarMolarVol[k];
}
}
if (TMolesInert > 0.0) {
if (m_gasPhase) {
double volI = TMolesInert * 8314.47215 * Temp / Pres;
m_totalVol += volI;
} else {
printf("unknown situation\n");
std::exit(-1);
}
m_totalVol = 0.0;
for (k = 0; k < NVolSpecies; k++) {
m_totalVol += PartialMolarVol[k] * Xmol[k];
}
m_totalVol *= v_totalMoles;
if (TMolesInert > 0.0) {
if (m_gasPhase) {
double volI = TMolesInert * 8314.47215 * Temp / Pres;
m_totalVol += volI;
} else {
printf("unknown situation\n");
std::exit(-1);
}
}
m_UpToDate_VolPM = true;
@ -895,18 +891,19 @@ namespace VCSnonideal {
/************************************************************************************/
/*
* updateLnActCoeffJac():
* _updateLnActCoeffJac():
*
*/
void vcs_VolPhase::updateLnActCoeffJac() {
void vcs_VolPhase::_updateLnActCoeffJac() {
int k, j;
double deltaMoles_j = 0.0;
/*
* Evaluate the current base activity coefficients.
*/
evaluateActCoeff();
* Evaluate the current base activity coefficients if necessary
*/
if (!m_UpToDate_AC) {
_updateActCoeff();
}
// Make copies of ActCoeff and Xmol for use in taking differences
std::vector<double> ActCoeff_Base(ActCoeff);
@ -939,7 +936,7 @@ namespace VCSnonideal {
* -> Note this calls setState_PX();
*/
_updateMoleFractionDependencies();
evaluateActCoeff();
_updateActCoeff();
/*
* Calculate the column of the matrix
*/
@ -962,7 +959,7 @@ namespace VCSnonideal {
*/
setMoleFractions(VCS_DATA_PTR(Xmol_Base));
_updateMoleFractionDependencies();
evaluateActCoeff();
_updateActCoeff();
}
/************************************************************************************/
@ -984,7 +981,8 @@ namespace VCSnonideal {
* mole number of species in the phase -> we always assume that
* they are out of date.
*/
updateLnActCoeffJac();
_updateLnActCoeffJac();
/*
* Now copy over the values
*/

View file

@ -146,15 +146,6 @@ namespace VCSnonideal {
void resize(int phaseNum, int numSpecies, const char *phaseName,
double molesInert = 0.0);
private:
//! Evaluate activity coefficients
/*!
* We carry out a calculation whenever UpTODate_AC is false. Specifically
* whenever a phase goes zero, we do not carry out calculations on it.
*/
void evaluateActCoeff() const;
public:
//! Evaluate activity coefficients and return the kspec coefficient
/*!
* We carry out a calculation whenever UpTODate_AC is false. Specifically
@ -239,19 +230,7 @@ namespace VCSnonideal {
* Units are potential
*/
double electricPotential() const;
private:
//! Gibbs free energy calculation for standard states
/*!
* Calculate the Gibbs free energies for the standard states
* The results are held internally within the object.
*
* @param TKelvin Current temperature
* @param pres Current pressure
*/
void GStar_calc() const;
public:
//! Gibbs free energy calculation for standard state of one species
/*!
* Calculate the Gibbs free energies for the standard state
@ -285,21 +264,6 @@ namespace VCSnonideal {
*/
double G0_calc_one(int kspec, double TKelvin);
private:
//! Molar volume calculation for standard states
/*!
* Calculate the molar volume for the standard states
* The results are held internally within the object.
*
* @param TKelvin Current temperature
* @param pres Current pressure
*
* Units are in m**3/kmol
*/
void VolStar_calc() const;
public:
//! Molar volume calculation for standard state of one species
/*!
* Calculate the molar volume for the standard states
@ -313,19 +277,8 @@ namespace VCSnonideal {
* @return molar volume of the kspec species's standard
* state (m**3/kmol)
*/
double VolStar_calc_one(int kglob, double TKelvin, double pres);
double VolStar_calc_one(int kglob) const;
private:
//! Calculate the partial molar volumes of all species and return the
//! total volume
/*!
* Calculates these quantitites internally
*
* @return total volume
*/
double VolPM_calc() const;
public:
//! Fill in the partial molar volume vector for VCS
/*!
* This routine will calculate the partial molar volumes for the
@ -362,23 +315,6 @@ namespace VCSnonideal {
*/
void setState_TP(double temperature_Kelvin, double pressure_PA);
private:
//! Evaluation of Activity Coefficient Jacobians
/*!
* This is the derivative of the ln of the activity coefficient
* with respect to mole number of jth species.
* (temp, pressure, and other mole numbers held constant)
*
* We employ a finite difference derivative approach here.
* Because we have to change the mole numbers, this is not
* a const function, even though the paradigm would say that
* it should be.
*
* @param moleNumbers Mole numbers are input.
*/
void updateLnActCoeffJac();
public:
// Downloads the ln ActCoeff jacobian into the VCS version of the
// ln ActCoeff jacobian.
/*
@ -438,7 +374,7 @@ namespace VCSnonideal {
* @param xmol Value of the mole fractions for the species
* in the phase. These are contiguous.
*/
void setMoleFractions (const double * const xmol);
void setMoleFractions(const double * const xmol);
//! Return a const reference to the mole fractions
const std::vector<double> & moleFractions() const;
@ -451,6 +387,59 @@ namespace VCSnonideal {
private:
//! Evaluate the activity coefficients at the current conditions
/*!
* We carry out a calculation whenever UpTODate_AC is false. Specifically
* whenever a phase goes zero, we do not carry out calculations on it.
*/
void _updateActCoeff() const;
//! Gibbs free energy calculation for standard states
/*!
* Calculate the Gibbs free energies for the standard states
* The results are held internally within the object.
*
* @param TKelvin Current temperature
* @param pres Current pressure
*/
void _updateGStar() const;
//! Molar volume calculation for standard states
/*!
* Calculate the molar volume for the standard states
* The results are held internally within the object.
*
* @param TKelvin Current temperature
* @param pres Current pressure
*
* Units are in m**3/kmol
*/
void _updateVolStar() const;
//! Calculate the partial molar volumes of all species and return the
//! total volume
/*!
* Calculates these quantitites internally
*
* @return total volume
*/
double _updateVolPM() const;
//! Evaluation of Activity Coefficient Jacobians
/*!
* This is the derivative of the ln of the activity coefficient
* with respect to mole number of jth species.
* (temp, pressure, and other mole numbers held constant)
*
* We employ a finite difference derivative approach here.
* Because we have to change the mole numbers, this is not
* a const function, even though the paradigm would say that
* it should be.
*
* @param moleNumbers Mole numbers are input.
*/
void _updateLnActCoeffJac();
//! Updates the mole fraction depenpencies
/*!
* Whenever the mole fractions change, this routine
@ -501,7 +490,7 @@ namespace VCSnonideal {
/*!
* The known types are listed at the top of this file.
*/
int EqnState;
int m_eqnState;
//! Number of element constraints within the problem
/*!

View file

@ -281,7 +281,7 @@ namespace VCSnonideal {
for (iphase = 0; iphase < NPhase; iphase++) {
Vphase = VPhaseList[iphase];
std::string EOS_cstr = string16_EOSType(Vphase->EqnState);
std::string EOS_cstr = string16_EOSType(Vphase->m_eqnState);
plogf("%16s %5d %5d %8d ", Vphase->PhaseName.c_str(),
Vphase->VP_ID, Vphase->SingleSpecies, Vphase->m_gasPhase);
plogf("%16s %8d %16e ", EOS_cstr.c_str(),

View file

@ -866,7 +866,7 @@ namespace VCSnonideal {
retn = VCS_PUB_BAD;
}
vPhase->EqnState = pub_phase_ptr->EqnState;
vPhase->m_eqnState = pub_phase_ptr->m_eqnState;
if (vPhase->NVolSpecies != pub_phase_ptr->NVolSpecies) {
plogf("%sNVolSpecies value have changed:%d %d\n", yo.c_str(),

View file

@ -230,7 +230,8 @@ VolStar_calc(int kglob, double TKelvin, double presPA)
if (UseCanteraCalls) {
AssertThrowVCS(m_VCS_UnitsFormat == VCS_UNITS_MKS, "Possible inconsistency");
int kspec = IndexSpeciesPhase;
vol = OwningPhase->VolStar_calc_one(kspec, TKelvin, presPA);
OwningPhase->setState_TP(TKelvin, presPA);
vol = OwningPhase->VolStar_calc_one(kspec);
} else {
switch(SSStar_Vol_Model) {
case VCS_SSVOL_CONSTANT: