20vcs_VolPhase::vcs_VolPhase(VCS_SOLVE* owningSolverObject) :
21 m_owningSolverObject(owningSolverObject)
25vcs_VolPhase::~vcs_VolPhase()
27 for (
size_t k = 0; k < m_numSpecies; k++) {
28 delete ListSpeciesPtr[k];
32void vcs_VolPhase::resize(
const size_t phaseNum,
const size_t nspecies,
33 const size_t numElem,
const char*
const phaseName,
34 const double molesInert)
36 AssertThrowMsg(nspecies > 0,
"vcs_VolPhase::resize",
"nspecies Error");
37 setTotalMolesInert(molesInert);
41 if (phaseNum == VP_ID_) {
42 if (strcmp(PhaseName.c_str(), phaseName)) {
44 "Strings are different: " + PhaseName +
" " +
45 phaseName +
" :unknown situation");
50 PhaseName = fmt::format(
"Phase_{}", VP_ID_);
52 PhaseName = phaseName;
56 m_singleSpecies =
false;
58 m_singleSpecies =
true;
61 if (m_numSpecies == nspecies && numElem == m_numElemConstraints) {
65 m_numSpecies = nspecies;
67 m_singleSpecies =
false;
70 IndSpecies.resize(nspecies,
npos);
72 if (ListSpeciesPtr.size() >= m_numSpecies) {
73 for (
size_t i = 0; i < m_numSpecies; i++) {
74 if (ListSpeciesPtr[i]) {
75 delete ListSpeciesPtr[i];
76 ListSpeciesPtr[i] = 0;
80 ListSpeciesPtr.resize(nspecies, 0);
81 for (
size_t i = 0; i < nspecies; i++) {
85 Xmol_.resize(nspecies, 0.0);
86 creationMoleNumbers_.resize(nspecies, 0.0);
87 creationGlobalRxnNumbers_.resize(nspecies,
npos);
88 for (
size_t i = 0; i < nspecies; i++) {
89 Xmol_[i] = 1.0/nspecies;
90 creationMoleNumbers_[i] = 1.0/nspecies;
91 if (IndSpecies[i] >= m_numElemConstraints) {
92 creationGlobalRxnNumbers_[i] = IndSpecies[i] - m_numElemConstraints;
94 creationGlobalRxnNumbers_[i] =
npos;
98 SS0ChemicalPotential.resize(nspecies, -1.0);
99 StarChemicalPotential.resize(nspecies, -1.0);
100 StarMolarVol.resize(nspecies, -1.0);
101 PartialMolarVol.resize(nspecies, -1.0);
102 ActCoeff.resize(nspecies, 1.0);
103 np_dLnActCoeffdMolNumber.resize(nspecies, nspecies, 0.0);
108 m_UpToDate_AC =
false;
109 m_UpToDate_VolStar =
false;
110 m_UpToDate_VolPM =
false;
111 m_UpToDate_GStar =
false;
112 m_UpToDate_G0 =
false;
117void vcs_VolPhase::elemResize(
const size_t numElemConstraints)
119 m_elementNames.resize(numElemConstraints);
120 m_elementActive.resize(numElemConstraints+1, 1);
122 m_formulaMatrix.resize(m_numSpecies, numElemConstraints, 0.0);
123 m_elementNames.resize(numElemConstraints,
"");
124 m_elemGlobalIndex.resize(numElemConstraints, npos);
125 m_numElemConstraints = numElemConstraints;
128void vcs_VolPhase::_updateActCoeff()
const
131 m_UpToDate_AC =
true;
134 TP_ptr->getActivityCoefficients(&ActCoeff[0]);
135 m_UpToDate_AC =
true;
138void vcs_VolPhase::_updateG0()
const
140 TP_ptr->getGibbs_ref(&SS0ChemicalPotential[0]);
141 m_UpToDate_G0 =
true;
144double vcs_VolPhase::G0_calc_one(
size_t kspec)
const
146 if (!m_UpToDate_G0) {
149 return SS0ChemicalPotential[kspec];
152void vcs_VolPhase::_updateGStar()
const
154 TP_ptr->getStandardChemPotentials(&StarChemicalPotential[0]);
155 m_UpToDate_GStar =
true;
158double vcs_VolPhase::GStar_calc_one(
size_t kspec)
const
160 if (!m_UpToDate_GStar) {
163 return StarChemicalPotential[kspec];
166void vcs_VolPhase::setMoleFractions(
const double*
const xmol)
169 for (
size_t k = 0; k < m_numSpecies; k++) {
173 if (std::fabs(sum) > 1.0E-13) {
174 for (
size_t k = 0; k < m_numSpecies; k++) {
178 _updateMoleFractionDependencies();
183void vcs_VolPhase::_updateMoleFractionDependencies()
186 TP_ptr->setMoleFractions(&Xmol_[m_MFStartIndex]);
187 TP_ptr->setPressure(Pres_);
189 if (!m_isIdealSoln) {
190 m_UpToDate_AC =
false;
191 m_UpToDate_VolPM =
false;
195const vector<double> & vcs_VolPhase::moleFractions()
const
200double vcs_VolPhase::moleFraction(
size_t k)
const
205void vcs_VolPhase::setMoleFractionsState(
const double totalMoles,
206 const double*
const moleFractions,
207 const int vcsStateStatus)
209 if (totalMoles != 0.0) {
213 throw CanteraError(
"vcs_VolPhase::setMolesFractionsState",
214 "inappropriate usage");
219 throw CanteraError(
"vcs_VolPhase::setMolesFractionsState",
220 "inappropriate usage");
225 m_vcsStateStatus = vcsStateStatus;
228 double fractotal = 1.0;
229 v_totalMoles = totalMoles;
230 if (m_totalMolesInert > 0.0) {
231 if (m_totalMolesInert > v_totalMoles) {
232 throw CanteraError(
"vcs_VolPhase::setMolesFractionsState",
233 "inerts greater than total: {} {}",
234 v_totalMoles, m_totalMolesInert);
236 fractotal = 1.0 - m_totalMolesInert/v_totalMoles;
239 for (
size_t k = 0; k < m_numSpecies; k++) {
240 Xmol_[k] = moleFractions[k];
241 sum += moleFractions[k];
244 throw CanteraError(
"vcs_VolPhase::setMolesFractionsState",
245 "inappropriate usage");
247 if (sum != fractotal) {
248 for (
size_t k = 0; k < m_numSpecies; k++) {
249 Xmol_[k] *= (fractotal /sum);
252 _updateMoleFractionDependencies();
255void vcs_VolPhase::setMolesFromVCS(
const int stateCalc,
256 const double* molesSpeciesVCS)
258 v_totalMoles = m_totalMolesInert;
260 if (molesSpeciesVCS == 0) {
261 AssertThrowMsg(m_owningSolverObject,
"vcs_VolPhase::setMolesFromVCS",
262 "shouldn't be here");
264 molesSpeciesVCS = &m_owningSolverObject->m_molNumSpecies_old[0];
266 molesSpeciesVCS = &m_owningSolverObject->m_molNumSpecies_new[0];
268 throw CanteraError(
"vcs_VolPhase::setMolesFromVCS",
"shouldn't be here");
270 }
else if (m_owningSolverObject) {
282 for (
size_t k = 0; k < m_numSpecies; k++) {
284 size_t kglob = IndSpecies[k];
285 v_totalMoles += std::max(0.0, molesSpeciesVCS[kglob]);
288 if (v_totalMoles > 0.0) {
289 for (
size_t k = 0; k < m_numSpecies; k++) {
291 size_t kglob = IndSpecies[k];
292 double tmp = std::max(0.0, molesSpeciesVCS[kglob]);
293 Xmol_[k] = tmp / v_totalMoles;
306 if (m_phiVarIndex !=
npos) {
307 size_t kglob = IndSpecies[m_phiVarIndex];
308 if (m_numSpecies == 1) {
309 Xmol_[m_phiVarIndex] = 1.0;
311 Xmol_[m_phiVarIndex] = 0.0;
313 double phi = molesSpeciesVCS[kglob];
314 setElectricPotential(phi);
315 if (m_numSpecies == 1) {
319 _updateMoleFractionDependencies();
320 if (m_totalMolesInert > 0.0) {
328 creationMoleNumbers_ = Xmol_;
333 m_vcsStateStatus = stateCalc;
336void vcs_VolPhase::setMolesFromVCSCheck(
const int vcsStateStatus,
337 const double* molesSpeciesVCS,
338 const double*
const TPhMoles)
340 setMolesFromVCS(vcsStateStatus, molesSpeciesVCS);
343 double Tcheck = TPhMoles[VP_ID_];
344 if (Tcheck != v_totalMoles) {
346 Tcheck = v_totalMoles;
348 throw CanteraError(
"vcs_VolPhase::setMolesFromVCSCheck",
349 "We have a consistency problem: {} {}", Tcheck, v_totalMoles);
354void vcs_VolPhase::updateFromVCS_MoleNumbers(
const int vcsStateStatus)
356 if ((!m_UpToDate || vcsStateStatus != m_vcsStateStatus) && m_owningSolverObject &&
358 setMolesFromVCS(vcsStateStatus);
362void vcs_VolPhase::sendToVCS_ActCoeff(
const int vcsStateStatus,
365 updateFromVCS_MoleNumbers(vcsStateStatus);
366 if (!m_UpToDate_AC) {
369 for (
size_t k = 0; k < m_numSpecies; k++) {
370 size_t kglob = IndSpecies[k];
371 AC[kglob] = ActCoeff[k];
375double vcs_VolPhase::sendToVCS_VolPM(
double*
const VolPM)
const
377 if (!m_UpToDate_VolPM) {
380 for (
size_t k = 0; k < m_numSpecies; k++) {
381 size_t kglob = IndSpecies[k];
382 VolPM[kglob] = PartialMolarVol[k];
387void vcs_VolPhase::sendToVCS_GStar(
double*
const gstar)
const
389 if (!m_UpToDate_GStar) {
392 for (
size_t k = 0; k < m_numSpecies; k++) {
393 size_t kglob = IndSpecies[k];
394 gstar[kglob] = StarChemicalPotential[k];
398void vcs_VolPhase::setElectricPotential(
const double phi)
401 TP_ptr->setElectricPotential(m_phi);
403 m_UpToDate_AC =
false;
404 m_UpToDate_VolStar =
false;
405 m_UpToDate_VolPM =
false;
406 m_UpToDate_GStar =
false;
409double vcs_VolPhase::electricPotential()
const
414void vcs_VolPhase::setState_TP(
const double temp,
const double pres)
416 if (Temp_ == temp && Pres_ == pres) {
419 TP_ptr->setElectricPotential(m_phi);
420 TP_ptr->setState_TP(temp, pres);
423 m_UpToDate_AC =
false;
424 m_UpToDate_VolStar =
false;
425 m_UpToDate_VolPM =
false;
426 m_UpToDate_GStar =
false;
427 m_UpToDate_G0 =
false;
430void vcs_VolPhase::setState_T(
const double temp)
432 setState_TP(temp, Pres_);
435void vcs_VolPhase::_updateVolStar()
const
437 TP_ptr->getStandardVolumes(&StarMolarVol[0]);
438 m_UpToDate_VolStar =
true;
441double vcs_VolPhase::VolStar_calc_one(
size_t kspec)
const
443 if (!m_UpToDate_VolStar) {
446 return StarMolarVol[kspec];
449double vcs_VolPhase::_updateVolPM()
const
451 TP_ptr->getPartialMolarVolumes(&PartialMolarVol[0]);
453 for (
size_t k = 0; k < m_numSpecies; k++) {
454 m_totalVol += PartialMolarVol[k] * Xmol_[k];
456 m_totalVol *= v_totalMoles;
458 if (m_totalMolesInert > 0.0) {
460 double volI = m_totalMolesInert *
GasConstant * Temp_ / Pres_;
463 throw CanteraError(
"vcs_VolPhase::_updateVolPM",
"unknown situation");
466 m_UpToDate_VolPM =
true;
470void vcs_VolPhase::_updateLnActCoeffJac()
472 double phaseTotalMoles = v_totalMoles;
473 if (phaseTotalMoles < 1.0E-14) {
474 phaseTotalMoles = 1.0;
478 if (!m_UpToDate_AC) {
484 TP_ptr->getdlnActCoeffdlnN(m_numSpecies, &np_dLnActCoeffdMolNumber(0,0));
485 for (
size_t j = 0; j < m_numSpecies; j++) {
486 double moles_j_base = phaseTotalMoles * Xmol_[j];
487 double*
const np_lnActCoeffCol = np_dLnActCoeffdMolNumber.ptrColumn(j);
488 if (moles_j_base < 1.0E-200) {
489 moles_j_base = 1.0E-7 * moles_j_base + 1.0E-13 * phaseTotalMoles + 1.0E-150;
491 for (
size_t k = 0; k < m_numSpecies; k++) {
492 np_lnActCoeffCol[k] = np_lnActCoeffCol[k] * phaseTotalMoles / moles_j_base;
496 double deltaMoles_j = 0.0;
498 vector<double> ActCoeff_Base(ActCoeff);
499 vector<double> Xmol_Base(Xmol_);
500 double TMoles_base = phaseTotalMoles;
503 for (
size_t j = 0; j < m_numSpecies; j++) {
506 double moles_j_base = phaseTotalMoles * Xmol_Base[j];
507 deltaMoles_j = 1.0E-7 * moles_j_base + 1.0E-13 * phaseTotalMoles + 1.0E-150;
511 phaseTotalMoles = TMoles_base + deltaMoles_j;
512 for (
size_t k = 0; k < m_numSpecies; k++) {
513 Xmol_[k] = Xmol_Base[k] * TMoles_base / phaseTotalMoles;
515 Xmol_[j] = (moles_j_base + deltaMoles_j) / phaseTotalMoles;
519 _updateMoleFractionDependencies();
523 v_totalMoles = TMoles_base;
530 setMoleFractions(&Xmol_Base[0]);
531 _updateMoleFractionDependencies();
535void vcs_VolPhase::sendToVCS_LnActCoeffJac(
Array2D& np_LnACJac_VCS)
539 _updateLnActCoeffJac();
542 for (
size_t j = 0; j < m_numSpecies; j++) {
543 size_t jglob = IndSpecies[j];
544 for (
size_t k = 0; k < m_numSpecies; k++) {
545 size_t kglob = IndSpecies[k];
546 np_LnACJac_VCS(kglob,jglob) = np_dLnActCoeffdMolNumber(k,j);
555 Pres_ = TP_ptr->pressure();
556 setState_TP(Temp_, Pres_);
557 m_phi = TP_ptr->electricPotential();
558 size_t nsp = TP_ptr->nSpecies();
559 size_t nelem = TP_ptr->nElements();
560 if (nsp != m_numSpecies) {
561 if (m_numSpecies != 0) {
562 warn_user(
"vcs_VolPhase::setPtrThermoPhase",
563 "Nsp != NVolSpeces: {} {}", nsp, m_numSpecies);
565 resize(VP_ID_, nsp, nelem, PhaseName.c_str());
567 TP_ptr->getMoleFractions(&Xmol_[0]);
568 creationMoleNumbers_ = Xmol_;
569 _updateMoleFractionDependencies();
573 m_isIdealSoln =
true;
575 m_isIdealSoln = TP_ptr->isIdeal();
584double vcs_VolPhase::totalMoles()
const
589double vcs_VolPhase::molefraction(
size_t k)
const
594void vcs_VolPhase::setCreationMoleNumbers(
const double*
const n_k,
595 const vector<size_t> &creationGlobalRxnNumbers)
597 creationMoleNumbers_.assign(n_k, n_k+m_numSpecies);
598 for (
size_t k = 0; k < m_numSpecies; k++) {
599 creationGlobalRxnNumbers_[k] = creationGlobalRxnNumbers[k];
603const vector<double>& vcs_VolPhase::creationMoleNumbers(
604 vector<size_t> &creationGlobalRxnNumbers)
const
606 creationGlobalRxnNumbers = creationGlobalRxnNumbers_;
607 return creationMoleNumbers_;
610void vcs_VolPhase::setTotalMoles(
const double totalMols)
612 v_totalMoles = totalMols;
613 if (m_totalMolesInert > 0.0) {
616 "vcs_VolPhase::setTotalMoles",
617 "totalMoles less than inert moles: {} {}",
618 totalMols, m_totalMolesInert);
620 if (m_singleSpecies && (m_phiVarIndex == 0)) {
623 if (totalMols > 0.0) {
632void vcs_VolPhase::setMolesOutOfDate(
int stateCalc)
635 if (stateCalc != -1) {
636 m_vcsStateStatus = stateCalc;
640void vcs_VolPhase::setMolesCurrent(
int stateCalc)
643 m_vcsStateStatus = stateCalc;
646bool vcs_VolPhase::isIdealSoln()
const
648 return m_isIdealSoln;
651size_t vcs_VolPhase::phiVarIndex()
const
653 return m_phiVarIndex;
656void vcs_VolPhase::setPhiVarIndex(
size_t phiVarIndex)
658 m_phiVarIndex = phiVarIndex;
660 if (m_singleSpecies && m_phiVarIndex == 0) {
667 return ListSpeciesPtr[kindex];
670int vcs_VolPhase::exists()
const
675void vcs_VolPhase::setExistence(
const int existence)
678 if (v_totalMoles != 0.0) {
680 "setting false existence for phase with moles");
682 }
else if (m_totalMolesInert == 0.0) {
683 if (v_totalMoles == 0.0 && (!m_singleSpecies || m_phiVarIndex != 0)) {
685 "setting true existence for phase with no moles");
690 "Trying to set existence of an electron phase to false");
692 m_existence = existence;
695size_t vcs_VolPhase::spGlobalIndexVCS(
const size_t spIndex)
const
697 return IndSpecies[spIndex];
700void vcs_VolPhase::setSpGlobalIndexVCS(
const size_t spIndex,
701 const size_t spGlobalIndex)
703 IndSpecies[spIndex] = spGlobalIndex;
704 if (spGlobalIndex >= m_numElemConstraints) {
705 creationGlobalRxnNumbers_[spIndex] = spGlobalIndex - m_numElemConstraints;
709void vcs_VolPhase::setTotalMolesInert(
const double tMolesInert)
711 if (m_totalMolesInert != tMolesInert) {
713 m_UpToDate_AC =
false;
714 m_UpToDate_VolStar =
false;
715 m_UpToDate_VolPM =
false;
716 m_UpToDate_GStar =
false;
717 m_UpToDate_G0 =
false;
718 v_totalMoles += (tMolesInert - m_totalMolesInert);
719 m_totalMolesInert = tMolesInert;
721 if (m_totalMolesInert > 0.0) {
723 }
else if (m_singleSpecies && (m_phiVarIndex == 0)) {
726 if (v_totalMoles > 0.0) {
734double vcs_VolPhase::totalMolesInert()
const
736 return m_totalMolesInert;
739size_t vcs_VolPhase::elemGlobalIndex(
const size_t e)
const
741 AssertThrow(e < m_numElemConstraints,
" vcs_VolPhase::elemGlobalIndex");
742 return m_elemGlobalIndex[e];
745void vcs_VolPhase::setElemGlobalIndex(
const size_t eLocal,
const size_t eGlobal)
748 "vcs_VolPhase::setElemGlobalIndex");
749 m_elemGlobalIndex[eLocal] = eGlobal;
752size_t vcs_VolPhase::nElemConstraints()
const
754 return m_numElemConstraints;
757string vcs_VolPhase::elementName(
const size_t e)
const
759 return m_elementNames[e];
765 for (
size_t k = 0; k < tPhase->
nSpecies(); k++) {
766 if (tPhase->
charge(k) != 0.0) {
786size_t vcs_VolPhase::transferElementsFM(
const ThermoPhase*
const tPhase)
796 ChargeNeutralityElement = ne;
803 if (ChargeNeutralityElement !=
npos) {
807 size_t eFound =
npos;
816 for (
size_t eT = 0; eT < nebase; eT++) {
819 m_elementActive[eT] = 0;
824 for (
size_t eT = 0; eT < nebase; eT++) {
831 if (eFound ==
npos) {
834 m_elementActive[ne] = 0;
836 m_elementNames[ne] = ename;
842 m_formulaMatrix.resize(ns, ne, 0.0);
847 for (
size_t eT = 0; eT < nebase; eT++) {
854 string pname = tPhase->
name();
856 pname = fmt::format(
"phase{}", VP_ID_);
858 e = ChargeNeutralityElement;
859 m_elementNames[e] =
"cn_" + pname;
862 for (
size_t k = 0; k < ns; k++) {
864 for (
size_t eT = 0; eT < nebase; eT++) {
865 m_formulaMatrix(k,e) = tPhase->
nAtoms(k, eT);
868 if (eFound !=
npos) {
869 m_formulaMatrix(k,eFound) = - tPhase->
charge(k);
874 for (
size_t k = 0; k < ns; k++) {
875 m_formulaMatrix(k,ChargeNeutralityElement) = tPhase->
charge(k);
882 if (ns == 1 && tPhase->
charge(0) != 0.0) {
890int vcs_VolPhase::elementType(
const size_t e)
const
892 return m_elementType[e];
895const Array2D& vcs_VolPhase::getFormulaMatrix()
const
897 return m_formulaMatrix;
900int vcs_VolPhase::speciesUnknownType(
const size_t k)
const
902 return m_speciesUnknownType[k];
905int vcs_VolPhase::elementActive(
const size_t e)
const
907 return m_elementActive[e];
910size_t vcs_VolPhase::nSpecies()
const
915string vcs_VolPhase::eos_name()
const
917 switch (m_eqnState) {
918 case VCS_EOS_CONSTANT:
920 case VCS_EOS_IDEAL_GAS:
922 case VCS_EOS_STOICH_SUB:
924 case VCS_EOS_IDEAL_SOLN:
926 case VCS_EOS_DEBEYE_HUCKEL:
927 return "Debeye Huckel";
928 case VCS_EOS_REDLICH_KWONG:
929 return "Redlich_Kwong";
930 case VCS_EOS_REGULAR_SOLN:
931 return "Regular Soln";
933 return fmt::format(
"UnkType: {:7d}", m_eqnState);
Header file for class ThermoPhase, the base class for phases with thermodynamic properties,...
A class for 2D arrays stored in column-major (Fortran-compatible) form.
Base class for exceptions thrown by Cantera classes.
size_t nSpecies() const
Returns the number of species in the phase.
double temperature() const
Temperature (K).
int elementType(size_t m) const
Return the element constraint type Possible types include:
double nAtoms(size_t k, size_t m) const
Number of atoms of element m in species k.
size_t nElements() const
Number of elements.
string elementName(size_t m) const
Name of the element with index m.
double charge(size_t k) const
Dimensionless electrical charge of a single molecule of species k The charge is normalized by the the...
string name() const
Return the name of the phase.
Base class for a phase with thermodynamic properties.
bool chargeNeutralityNecessary() const
Returns the chargeNeutralityNecessity boolean.
Properties of a single species.
#define AssertThrow(expr, procedure)
Assertion must be true or an error is thrown.
#define AssertThrowMsg(expr, procedure,...)
Assertion must be true or an error is thrown.
const double GasConstant
Universal Gas Constant [J/kmol/K].
void warn_user(const string &method, const string &msg, const Args &... args)
Print a user warning raised from method as CanteraWarning.
Namespace for the Cantera kernel.
const size_t npos
index returned by functions to indicate "no position"
static bool chargeNeutralityElement(const ThermoPhase *const tPhase)
This utility routine decides whether a Cantera ThermoPhase needs a constraint equation representing t...
static bool hasChargedSpecies(const ThermoPhase *const tPhase)
This function decides whether a phase has charged species or not.
bool vcs_doubleEqual(double d1, double d2)
Simple routine to check whether two doubles are equal up to roundoff error.
Contains declarations for string manipulation functions within Cantera.
Header for the object representing each phase within vcs.
#define VCS_SPECIES_TYPE_INTERFACIALVOLTAGE
Unknown refers to the voltage level of a phase.
#define VCS_STATECALC_OLD
State Calculation based on the old or base mole numbers.
#define VCS_STATECALC_TMP
State Calculation based on a temporary set of mole numbers.
#define VCS_PHASE_EXIST_NO
Phase doesn't currently exist in the mixture.
#define VCS_PHASE_EXIST_ALWAYS
#define VCS_ELEM_TYPE_ABSPOS
Normal element constraint consisting of positive coefficients for the formula matrix.
#define VCS_STATECALC_NEW
State Calculation based on the new or tentative mole numbers.
#define VCS_SPECIES_TYPE_MOLNUM
Unknown refers to mole number of a single species.
#define VCS_ELEM_TYPE_ELECTRONCHARGE
This refers to conservation of electrons.
#define VCS_ELEM_TYPE_CHARGENEUTRALITY
This refers to a charge neutrality of a single phase.
#define VCS_PHASE_EXIST_YES
Phase is a normal phase that currently exists.
#define VCS_PHASE_EXIST_ZEROEDPHASE
Phase currently is zeroed due to a programmatic issue.
Header file for the internal object that holds the vcs equilibrium problem (see Class VCS_SOLVE and C...